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1.  Summary 


We  develop  a  membrane  surface  wave  code  named  SurfMembrane  for  the  purpose  of  predicting 
surface  wave  amplitudes  in  heterogeneous  structures.  The  code  models  diffraction  and  scattering 
in  heterogeneous  structures  and  runs  orders  of  magnitude  faster  than  a  corresponding  3D  finite 
difference  calculation.  The  main  limitation  of  the  membrane  calculation  is  that  it  is  restricted  to  a 
narrow  frequency  band,  small  enough  that  the  frequency  dependence  of  phase  velocity  can  be 
ignored.  We  compared  the  membrane  calculations  with  a  finite  difference  calculation  for  the  same 
complex  structure  and  found  good  agreement  in  amplitude  variations. 

We  ran  SurfMembrane  for  73  earthquakes  recorded  on  the  US  Array  using  the  CMT  solutions  for 
those  events  to  generate  the  source,  and  compared  the  calculations  with  data.  We  found  that  while 
the  membrane  calculations  produced  amplitude  variations  comparable  to  those  observed,  the 
detailed  results  are  very  sensitive  to  small  errors  in  both  the  earth  model  and  source  radiation 
pattern,  so  it  was  not  possible  to  generate  predicted  amplitude  changes  sufficiently  accurate  to 
correct  amplitudes  at  a  specific  location. 

We  added  to  SurfMembrane  the  capability  to  calculate  sensitivity  kernels  which  predict  variations 
in  travel  time  as  a  function  of  earth  structure,  and  we  compared  the  results  from  the  sensitivity 
kernels  to  ray  theory  predictions.  The  sensitivity  kernels  can  predict  changes  due  to  off-axis 
structural  variations  that  cannot  be  modeled  by  ray  theory,  however  for  most  cases  we  found  that 
predictions  from  ray  theory  along  great  circle  paths  gave  more  accurate  predictions  than  the 
sensitivity  kernels.  This  is  most  likely  due  to  the  Bom  approximation  used  in  calculation  of 
sensitivity  kernels,  which  is  only  valid  for  small  perturbations. 

We  inverted  a  large  Eurasian  and  northern  African  data  set  for  Q  structure.  The  data  set  retrieved 
from  IRIS  consisted  of  59,000  waveforms  from  1850  Eurasian  and  African  earthquakes,  all  with 
CMT  solutions.  These  were  added  to  the  smaller  data  set  used  in  the  Stevens  et  al  (2008)  study. 
Both  the  new  and  the  old  data  sets  were  subjected  to  an  intensive  quality  control  procedure  to 
remove  bad  or  questionable  data,  deep  events  and  events  with  inaccurate  CMT  solutions.  The  final 
data  set  contained  23,148  waveforms  from  998  events.  Inversion  results  are  similar  to  the  2008 
study  and  cover  a  larger  region  of  Eurasia.  There  is  a  band  of  high  attenuation  stretching  across 
the  Middle  East  from  the  Mediterranean  Sea  to  India. 

We  ran  SurfMembrane  calculations  for  88  well-distributed  earthquakes  from  the  Eurasian  data  set 
and  compared  the  results  with  data.  We  find  that  the  SurfMembrane  calculations  can  be  used  to 
correct  Ms  for  diffraction  out  to  distances  of  about  2000  km,  however  at  greater  distances 
diffraction  is  too  sensitive  to  small  variations  in  earth  structure  to  allow  reliable  predictions. 

Together  with  this  report,  we  are  delivering  SurfMembrane  3.0,  together  with  test  cases  that  can 
be  run  to  ensure  the  program  is  working  correctly.  The  program  manual,  which  includes 
instructions  for  running  the  program  and  descriptions  of  all  file  formats,  is  given  in  the  appendices 
to  this  report.  We  are  also  delivering  the  earth  structures  and  derived  attenuation  coefficients  from 
the  Eurasian  Q  inversion  study. 
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2.  Introduction 


The  objective  of  this  project  is  to  improve  modeling  and  prediction  of  surface  wave  amplitudes  at 
short  periods  (8-20  seconds)  for  nuclear  monitoring  purposes.  This  follows  previous  projects 
which  developed  global  surface  wave  dispersion  maps  on  a  one  degree  grid  (Stevens  et  al,  2005), 
and  inverted  surface  wave  data  in  Eurasia  for  attenuation  (Stevens  et  al,  2008).  The  second  project 
showed  that  attenuation  inversions  were  strongly  affected  by  amplitude  variations  caused  by  earth 
structure,  and  so  in  this  project  we  are  implementing  and  testing  a  new  procedure  based  on 
membrane  surface  waves  to  correct  for  structural  variations  in  amplitude,  and  perform 
comparisons  with  both  finite  difference  calculations  and  data  to  see  if  these  structural  effects  can 
be  predicted  and  removed. 

We  developed  a  membrane  surface  wave  code  named  “SurfMembrane.”  As  a  first  test, 
SurfMembrane  calculations  were  compared  with  the  results  of  a  3D  finite  difference  calculation 
corresponding  to  the  complex  region  between  Lop  Nor  and  KNET  (Stevens  et  al,  2008).  Phase 
velocities  for  several  frequencies  were  derived  from  the  KNET  model  of  velocity  and  density  vs. 
depth  at  each  point  in  the  3D  grid.  SurfMembrane  was  run  for  each  frequency  and  compared  with 
results  from  the  3D  calculation  filtered  in  approximately  the  same  frequency  band.  Amplitude 
corrections  derived  from  the  SurfMembrane  calculations  significantly  improved  the  consistency 
of  the  results.  This  shows  that  if  the  source  and  structure  are  well-defined,  SurfMembrane  can  be 
used  to  generate  structural  amplitude  corrections. 

SurfMembrane  was  then  applied  to  US  Array  data  and  Eurasian  surface  wave  data,  using  the  best 
available  phase  velocity  models  for  each  and  calculating  surface  waves  for  events  with  known 
CMT  solutions.  As  discussed  in  the  following  sections,  the  SurfMembrane  calculations  generally 
predict  the  magnitude  of  the  expected  variability,  but  there  is  too  much  sensitivity  to  the  details  of 
the  earth  structure  to  make  point  by  point  predictions.  For  earthquakes,  we  also  find  that  small 
errors  in  radiation  pattern  derived  from  CMT  solutions  cause  amplitude  variations  that  dominate 
over  the  diffraction  effects. 

We  also  use  SurfMembrane  to  calculate  sensitivity  kernels,  which  predict  the  change  in  travel  time 
as  a  function  of  variations  in  earth  structure,  including  variations  on  earth  structure  that  do  not  lie 
along  the  great  circle  path.  We  perform  a  detailed  comparison  of  the  predictions  from  the 
sensitivity  kernels  with  the  predictions  from  ray  theory  along  great  circle  paths. 

We  extend  the  inversion  of  surface  wave  attenuation  for  Q  structure  that  was  performed  in  2008, 
adding  much  more  data  and  extending  the  range  into  northern  Africa,  we  retrieved  59,000 
waveforms  from  1850  Eurasian  and  African  events.  We  also  retrieved  the  locations,  depths  and 
origin  times  for  all  of  these  events  from  the  International  Seismological  Center  (ISC),  and  replaced 
the  CMT  origins  and  depths  with  the  better  ISC  origins.  These  data  were  then  combined  with  the 
data  set  from  the  2008  study  and  put  through  an  aggressive  data  quality  procedure  to  prune  out  bad 
or  noisy  data,  and  deep  events  or  data  with  bad  depths  or  focal  mechanisms.  This  resulted  in  a  high 
quality  data  set  of  23,148  spectra  from  998  events.  There  were  an  average  of  14.8  spectral  data 
points  for  each  waveform,  leading  to  a  total  data  set  of  343,322  data  points. 

Appendix  A  is  a  User  Manual  for  the  SurfMembrane  code,  and  Appendix  B  and  C  are  a  description 
of  the  code  and  data  delivered  together  with  this  report. 
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3.  Technical  Approach 
3.1  Surface  Wave  Data  Set 

We  assembled  two  large  data  sets  of  surface  waves,  all  from  events  large  enough  to  have  CMT 
solutions: 

1.  Areas  in  the  continental  United  States  where  most  of  the  source  to  receiver  path  is 
covered  by  the  US  Array.  There  were  73  such  events,  each  recorded  by  up  to  500  US 
Array  stations  (Figure  1). 


45  N 


30  N 


120°  W  90°  W  60°  W 

Figure  1.  Events  (triangles)  that  occurred  while  the  US  Array  was  passing  over  them  and  were 
large  enough  to  have  Centroid  Moment  Tensor  (CMT)  solutions  are  shown  as  solid  triangles.  US 
Array  stations  are  shown  as  +  symbols,  colored  to  show  stations  that  were  operating  at  the  same  time. 


2.  Eurasian  and  North  African  earthquakes  recorded  on  GSN  stations.  We  collected  59,000 
waveforms  from  1850  Eurasian  and  African  events  (Figure  2).  Measurements  made  from 
this  data  set  were  added  to  the  data  set  used  in  the  previous  project,  and  an  extensive 
quality  control  procedure  was  applied  to  the  entire  data  set  (see  section  4.6). 
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Figure  2.  Eurasian  and  African  GSN  data  used  in  this  project.  Top:  Red  events  are  in  the  region  of 
interest  with  predominately  continental  paths;  blue  events  are  outside  the  region  of  interest.  Bottom: 
GSN  stations  that  recorded  these  events.  Yellow  stations  are  within  the  region  of  interest. 
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3.2  Measurement  of  Surface  Wave  Amplitudes 

All  new  data  for  this  project  was  retrieved  from  the  Incorporated  Research  Institutions  for 

Seismology  (IRIS)  using  the  Matlab  irisFetch  routine.  Two  types  of  measurements  were  made: 

1 .  Data  was  band-pass  filtered  with  a  3  pole  two-pass  Butterworth  filter  over  a  frequency  band 
of/ ±0.1/,  for  a  range  of  frequencies ,  /  corresponding  to  periods  from  8  to  25  seconds. 
Amplitudes  were  measured  from  the  peak  of  the  envelope  function  derived  from  the  Hilbert 
transform  of  the  waveform.  Amplitudes  were  corrected  for  instrument  response  using  the 
instrument  responses  also  retrieved  from  IRIS. 

2.  A  phase-matched  filter  was  constructed  using  a  synthetic  seismogram  calculated  using  the 
CMT  solution  and  structure  and  attenuation  models  from  Stevens  et  al  (2008).  Phase  match 
filtered  seismograms  were  windowed  at  ±100  seconds,  Fourier  transformed  and  corrected 
for  instrument  response.  Smoothed  spectra  were  then  used  for  Q  inversion  (section  4.5). 

3.3  Earth  Models 


We  are  using  two  sets  of  earth  models:  1)  Leidos  global  models  developed  using  surface  waves  on 
a  1  degree  grid  (Stevens  et  al,  2005);  and  2)  a  higher  resolution  model  (1/4  degree)  for  the  Western 
United  States  (WUS)  from  the  University  of  Colorado  developed  from  ambient  noise  and 
earthquake  data  (Shen  et  al,  2013).  Phase  velocity  maps  at  10  seconds  over  the  Western  United 
States  are  shown  in  Figure  3. 
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Figure  3.  Left:  Leidos  1-degree  phase  velocity  model  of  the  Western  United  States  at  10  seconds 
period.  Right:  CU  ^-degree  phase  velocity  model  of  the  same  area  at  10  seconds.  Distances  are  in 
km,  phase  velocity  in  m/s.  Point  (0,0)  corresponds  to  an  event  near  Reno,  NV  at  Lat  39.520,  Lon  - 
119.930. 
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3.4  Membrane  Surface  Wave  Method 


The  membrane  surface  wave  technique  was  first  developed  by  Tanimoto  (1990),  and  used  to 
model  the  complex  propagation  of  surface  waves  from  the  Whittier  earthquake  (Figure  4).  The 
membrane  surface  wave  technique  reduces  a  3 -dimensional  elastic  problem  to  a  two- 
dimensional  scalar  wave  problem,  with  the  surface  wave  modeled  as  a  wave  traveling  through  a 
surface  with  variable  phase  velocity.  Since  phase-velocity  is  frequency  dependent,  the 
technique  is  only  expected  to  be  valid  over  a  narrow  enough  frequency  band  that  the  phase 
velocity  does  not  change  significantly.  The  advantages  of  the  membrane  surface  wave 
technique  are  that  1)  it  is  orders  of  magnitude  faster  than  3-dimensional  finite-difference;  and 
2)  it  accurately  handles  diffraction  and  scattering  along  a  complex  travel  path. 


Figure  4.  Membrane  surface  wave  calculation 
from  Tanimoto  (1990). 
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4.  Results  and  Discussion 


4.1  Implementation  of  SurfMembrane 

We  developed  a  membrane  surface  wave  code  named  “SurfMembrane”  using  a  finite  difference 
equation  solver  (Tape,  2003)  and  absorbing  boundaries  (Clayton  and  Engquist,  1977).  The 
membrane  surface  wave  technique  can  be  applied  to  either  a  planar  or  spherical  surface.  Here,  we 
have  applied  it  to  a  planar  surface.  In  either  case,  we  need  to  solve  the  wave  equation: 


1  d2u 


-  V2u  =  / 


(1) 


where  phase  velocity  c  =  c(x,  y )  ,  wave  field  u  =  u(x,  y,t) ,  forcing  function  /  =  /( x,  y,t ) ,  and 
82u  d2u 

V2w  =  — -  +  — -  is  the  Laplacian  of  u.  The  second-order  leap-frog  finite-difference  approximation 

dx  dy 

of  the  wave  equation  given  explicitly  in  Tape  (2003)  is  remarkably  simple  and  straightforward  to 
implement.  The  Clayton  and  Enquist  (1977)  absorbing  boundaries  do  an  excellent  job  of  removing 
any  reflected  waves  from  the  boundaries. 


The  wave  field  can  be  initiated  with  either  a  forcing  function,  an  initial  wave  field  or  by  specifying 
the  wave  field  on  a  closed  surface  (in  the  latter  case,  the  wave  field  is  only  valid  outside  the  closed 
surface.  Since  we  want  to  simulate  the  motion  of  a  surface  wave  from  an  explosion  or  earthquake, 
we  developed  an  initialization  procedure  that  accurately  reproduces  the  surface  wave  by 
calculating  a  fundamental  mode  surface  wave  in  the  structure  close  to  the  source  at  a  number  of 
points  a  few  km  away  from  and  encircling  the  source  point.  We  then  drive  the  membrane  solver 
with  the  motion  defined  at  these  points,  band-pass  filtered  so  that  the  frequency  corresponds 
closely  to  the  frequency  used  in  the  phase  velocity  map.  To  model  a  particular  event  with  a  CMT 
solution,  the  surface  wave  is  generated  from  the  CMT  solution  for  the  event. 

An  example  is  shown  in  Figure  5  for  a  strike  slip  earthquake  at  a  depth  of  10  km.  Surface  wave 
motion  was  calculated  using  a  modal  code  at  40  points  on  a  square  20  km  across  centered  on  the 
source  and  then  narrow-band  filtered  near  10  seconds  period.  SurfMembrane  successfully 
propagated  that  motion  and  preserved  the  initial  radiation  pattern. 


Figure  5.  Left:  Initial  motion  from  a  strike-slip  earthquake  on  a  closed  square  around  the  source. 
Right:  motion  at  170  seconds  showing  that  the  radiation  pattern  is  preserved. 
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We  found  that  it  is  necessary  to  calculate  the  exact  values  of  the  Hankel  functions  in  generating 
the  surface  wave  drive  functions  (See  Harkrider  et  al,  1994,  equation  51).  Use  of  the  asymptotic 
form  of  the  Hankel  functions,  which  is  common  in  surface  wave  codes,  introduces  an  inaccurate 
source  phase  that  leads  to  some  very  odd  radiation  patterns,  particularly  when  there  are  both  dip- 
slip  and  strike-slip  components  to  the  source  (See  Figure  6). 
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Figure  6.  Left:  Radiation  pattern  of  event  with  strike  slip  and  dip  slip  components  using  the  full 
Hankel  function  equations  for  the  source.  Right:  same  using  the  asymptotic  approximation.  Both 
calculations  are  at  20  seconds. 

4.2  Test  of  SurfMembrane  on  Embedded  Tarim  Basin  Structure 


SurfMembrane  was  tested  by  comparison  with  a  3D  finite  difference  calculation  of  propagation  of 
explosion-generated  wave  motion  across  the  Tarim  Basin  (Stevens  et  al,  2008).  The  earth  model 
has  a  rectangular  region  with  the  Tarim  Basin  structure  (Table  2)  embedded  in  a  Eurasian  Shield 
structure  (Table  1).  Both  models  are  derived  from  the  global  earth  models  of  Stevens  et  al  (2005). 


Table  1.  Asian  Shield  earth  model 


Depth  (km) 

Vp  (km/s) 

Vs  (km/s) 

Density  (g/cm3) 

10.5 

6.148 

3.451 

2.643 

20.5 

6.224 

3.494 

2.671 

50.0 

6.800 

3.817 

2.881 

66.7 

8.042 

4.514 

3.334 

83.3 

8.042 

4.514 

3.334 

100 

8.247 

4.629 

3.409 

00 

8.404 

4.717 

3.466 

Table  2.  Tarim  Basin  earth  model 

Depth  (km) 

Vp  (km/s) 

Vs  (km/s) 

Density  (g/cm3) 

6.5 

4.600 

2.550 

2.480 

10.0 

5.396 

3.029 

2.369 

20.0 

5.496 

3.085 

2.405 

30.0 

6.428 

3.608 

2.745 

40.0 

6.802 

3.818 

2.882 

50.0 

7.349 

4.125 

3.081 

68.7 

7.576 

4.253 

3.164 

83.3 

8.042 

4.514 

3.334 

100 

8.247 

4.629 

3.409 

00 

8.404 

4.717 

3.466 
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While  the  finite  difference  calculation  uses  the  full  3D  structures,  the  membrane  surface  wave 
calculation  uses  only  the  phase  velocity  at  a  particular  period.  At  10  seconds  period,  the  phase 
velocities  are  3.264  km/s  and  2.688  km/s  for  the  Eurasian  and  Tarim  Basin  structures,  respectively. 
A  snapshot  of  the  vertical  velocity  at  240  seconds  is  shown  in  Figure  7.  The  results  are  very  similar, 
showing  strong  diffraction  around  the  inclusion  in  both  cases. 


at  240  seconds  x  Time  240, 1 497,  Max  0. 144 


Figure  7.  Left:  3D  finite  difference  calculation  of  vertical  velocity  from  an  explosion  propagating 
past  a  low  velocity  inclusion  (marked  with  square)  modeled  after  the  Tarim  Basin  (from  Stevens  et 
al,  2008).  Right:  Same  calculation  with  membrane  surface  wave  code  for  10  second  period.  Axes 
show  distance  east  and  north  in  km.  The  explosion  source  is  at  North  0 \  East  1150. 
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4.3  Calculation  and  Evaluation  of  Membrane  Surface  Wave  Finite  Frequency 
Sensitivity  Kernels 

Finite-frequency  sensitivity  kernels  (e.g.  Tromp  et  al.  2005,  Chen  et  al  2007)  give  an  estimate  of 
the  change  in  an  observable,  such  as  arrival  time,  due  to  a  change  in  an  earth  model.  A  great  deal 
of  research  has  been  performed  over  the  past  two  decades  on  the  construction  of  such  kernels  and 
their  use  in  inversion  for  2D  and  3D  earth  structure.  In  3D  elastic  structures,  these  kernels  can  be 
quite  expensive  and  time  consuming  to  calculate,  but  for  membrane  surface  waves  they  can  be 
calculated  very  quickly. 

Our  interests  here  are  to  1)  show  how  membrane  sensitivity  kernels  can  be  calculated  using  the 
SurfMembrane  code,  and  2)  assess  their  usefulness  relative  to  traditional  ray  theory  tomography 
for  performing  inversions  for  earth  structure.  In  ray  theory  tomography,  travel  times  are 
determined  by  integrating  along  a  great  circle  (or  modified  great  circle)  path,  while  in  finite- 
frequency  tomography,  travel  times  are  determined  by  integrating  over  the  earth’s  surface. 

4.3.1  Calculation  of  sensitivity  kernels 

The  derivation  below  closely  follows  Peter  et  al  (2007),  and  also  follows  Tape  et  al  (2007)  and 
Tromp  et  al  (2005).  The  membrane  surface  wave  technique  solves  the  wave  equation: 


1  d2u 

-flit2 


-V2u  =  f 


(2) 


where  phase  velocity  c  =  c(x,  y )  ,  wave  field  u  =  u(x,  y,t )  and  forcing  function  /  =  fix,  y,t )  and 


V2u 


d2u  d2u 
dx2  dy2 


is  the  Laplacian  of  u.  The  Green’s  function  is  defined  as  the  solution  to: 


c2  dt 2 


-V2 


G(x,x  =  S(x  -x 


(3) 


The  variation  in  the  waveform  u  due  to  a  perturbation  in  the  phase  velocity  c  is  given  by: 

2 


Su(x,t)=  [  [  r-—G(x,x';t-t') 
nc  (*  ) 


d2u(x',t')  Sc(x ') 

K  ’  V  ’  dx'dt' 


dr 


The  variation  in  travel  time  to  a  receiver  at  xr  is  given  by: 

1  cT  ,.du(xr,t) 


ST  = — [  w(t ) — -  —  -Su(x  ,t)dt 
N  Jo  dt  X  r  ’ 


Where  w(t)  is  a  windowing  function  and  the  receiver  dependent  normalization  N  is: 


N, 


c1  { 

=  Jo 


dt2 


Note  that  N.  will  always  be  negative.  Combining  (4)  with  (5): 


(4) 


(5) 


(6) 
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1  T 

ST  = — f  wit) 

N  Jo 


du(xr,t) 


which  can  be  written: 


dt 


,  d2u(x  ',t ')  Sc(x ') 


dt2 


dx'  dt'  dt 


(7) 


c  1  pt  2  d2u(x \t')  Sc(x ')  pT-t  ,  ,  ,  . 

8T=  f — f  - V  /  ;  /  /  G{x\xr\T-t-t')w{T-t') 

{  nJ °  c2(x')  at2  c(x')  Jo  v  r  ;  v  ; 


duix.T -t') 

— - ’-dx'dt'dt 


dt 


(8) 


The  last  integral  corresponds  to  a  field  generated  by  a  disturbance  at  the  receiver  point,  so  we 
define  this  “adjoint  field”  as: 


1  r _ _  .  . t\  t m  .a du{xr,T-t') 


id  {x,xr\T  -t)  =  —~  £  G[x,xr;T -t-t')w(T -t')- 


dt 


dt' 


Which  reduces  equation  8  to: 


r  tT  2  d2u[x t)  Sc(x ')  +/  , 

<57’=ff  -^7- — - \  ;  /  Ad(T,xr\T-t)dx'dt 

c2(x')  dt2  c(x')  V  r  ’ 


Where  id  is  the  adjoint  field  generated  by  the  “adjoint  source”: 


z1 


yjN  dt  v  rJ 


We  define  the  sensitivity  kernel  as: 


d2u(yx,t ) 


dt2 


dt 


And  so: 


ST  =\^j^-K(x,xr)dx 

a  c\x ) 


(9) 


(10) 


(11) 


(12) 


(13) 


Note  that  we  have  chosen  the  sign  convention  that  ST  =  Tobs  -Tcalc  so  the  time  residual  is  positive 
for  a  later  than  predicted  arrival  time. 


There  is  one  remaining  problem:  delta  functions  do  not  exist  in  the  real  world  and  so  we  cannot 
implement  equation  1 1  directly.  Instead,  we  integrate  equation  1 1  over  the  grid  cell  containing  the 
receiver  with  area  Ar  and  define  the  alternate  quantities  (Peter  et  al.  refer  to  these  as  “discretized” 
quantities): 


— +  .  .  1  du(x,T-t) 

f(x,t)  =  —  w(T-t)—d^ - 

y  ’  N„  dt 


(14) 


M+(x,xr;r-f) 


—  [T  ’ G(x,x  ;T -t- 
Nr  Jo  V 


t')w(T-t') 


du(xr,T 

dT 


(15) 
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Then: 


, _ \  2  rr_+/ N d2u(x,t ), 

K(X’X’^c2(x)aJ«U  < x’x''T-‘ )  d‘ 


A  similar  analysis  from  Tape  et  al  (2007)  expresses  the  kernel  in  terms  of  spatial  gradients: 

r--x  -2 


K(x,xr)  =  —  £  Vwf  (x,x/,T  -t)(ftu(x,tyit 


(16) 


(17) 


Equations  16  and  17  give  identical  results. 

m"1'  is  the  quantity  calculated  by  the  membrane  code  in  response  to  the  force  in  equation  14,  and 
the  sensitivity  kernel  is  then  calculated  from  equation  16.  The  adjoint  force  (14)  is  derived  from 
the  time  reversed  displacement  at  the  receiver,  which  is  then  reversed  in  time  again  to  calculate 
the  kernel  (16). 

SurfMembrane  has  been  modified  to  do  this  entire  calculation.  To  calculate  the  kernel  using 
SurfMembrane,  run  the  forward  calculation  saving  the  displacement  at  receiver  points,  set  the  flag 
“kernel_from_stadata”  and  the  output  file  “kernelfilename”.  SurfMembrane  first  runs  through  the 
forward  calculation  and  calculates  the  displacement  at  receiver  points.  It  then  differentiates  and 
normalizes  the  receiver  displacements  and  reverses  the  resulting  velocity  in  time  (equation  14).  It 
then  runs  a  forward  calculation  for  each  receiver  point,  calculating  the  displacement  field 
u f  (x,  xr;t) ,  which  is  then  numerically  integrated  with  the  initial  acceleration  field  (read  back  from 

the  initial  calculation  and  double  differentiated)  to  calculate  the  kernel  using  equation  16.  Note 
that  the  entire  field  for  all  times  must  be  saved  for  the  initial  source,  but  it  is  not  necessary  to  save 
the  calculation  for  the  adjoint  source  as  this  is  calculated  and  integrated  on  the  fly. 

4.3.2  Comparison  of  sensitivity  kernels  with  ray  theory 

For  comparison,  ray  theory  only  considers  variations  in  velocity  along  the  ray  path  between  the 
source  and  receiver.  In  that  case  the  variation  in  travel  time  (for  small  perturbations)  is  given  by: 


AT  =  - 


A  r  Ac 
c  c 


(18) 


where  A r  is  the  linear  distance  over  which  the  perturbation  occurs  and  Ac/c  is  the  fractional 
perturbation  in  phase  velocity.  This  compares  with  equation  13,  in  which  the  variation  in  travel 
time  is  related  to  a  surface  integral  with  the  kernel. 

Figure  8  shows  the  10  second  period  sensitivity  kernel  calculated  for  a  structure  with  a  uniform 
velocity  of  3264  m/s,  between  a  source  and  receiver  separated  by  1600  km.  Note  that  although  the 
great  circle  path  between  the  source  and  receiver  shows  an  increase  in  travel  time  with  a  decrease 
in  velocity,  consistent  with  ray  theory,  a  perturbation  120  km  off  axis  produces  a  change  of  the 
opposite  sign,  so  a  low  velocity  zone  at  this  location  is  predicted  to  cause  a  decrease  in  travel  time. 
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Figure  8.  Left:  Sensitivity  kernel  calculated  between  the  source  (black  circle)  and  station  (yellow 
circle)  in  a  uniform  structure.  The  kernel  is  calculated  at  a  period  of  10  seconds  with  units  of  m 2 
corresponding  to  the  change  in  fractional  travel  time  per  change  in  fractional  phase  velocity  per 
square  meter.  Grid  units  are  km.  Right:  Cross  section  across  the  center  of  the  kernel. 

To  compare  with  ray  theory  directly,  we  calculated  the  change  in  travel  time  for  vertical  strips 
across  this  kernel,  calculating  equation  13  with  the  kernel  above  and  equation  18  for  ray  theory. 
The  strips  used  and  the  resulting  travel  time  changes  are  shown  in  Figure  9.  The  results  for  the 
kernel  area  integral  and  ray  line  integral  are  almost  identical  in  this  case,  except  that  the  ray  theory 
time  delay  is  1.7%  longer. 
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Figure  9.  Comparison  of  travel  times  calculated  by  ray  theory  and  kernel  integral  across  each  patch 
shown  on  the  left.  Time  delay  units  are  seconds  per  phase  velocity  fraction,  so  a  1%  change  in  phase 
velocity  in  a  patch  results  in  1%  of  the  time  delay  shown  above. 

To  further  check  the  accuracy  of  these  predictions,  we  ran  three  test  cases  with  perturbations  at  the 
locations  shown  in  Figure  10.  In  each  case  the  perturbation  was  a  region  with  a  low  velocity  of 
2688  m/s  (-17.6%  perturbation).  The  three  cases  correspond  to  1)  a  small  region  along  the  direct 
path  between  the  source  and  receiver;  2)  a  small  region  in  the  area  expected  to  produce  a  travel 
time  change  of  the  opposite  sign  from  the  direct  path;  and  3)  a  region  completely  crossing  the 
kernel. 
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Figure  10.  Calculations  were  performed  with  low  velocity  perturbations  in  the  3  regions  shown 
above. 

To  calculate  the  change  in  travel  time,  the  waveforms  from  the  perturbed  cases  were  cross- 
correlated  with  the  waveform  from  the  uniform  structure  and  then  Hilbert  transformed  to  find  the 
peak  arrival  time  (Figure  11).  Note  that  both  of  these  perturbations  produce  a  strong  diffracted 
wave  that  interferes  with  the  direct  wave  and  produces  a  secondary  arrival.  Travel  times  are  derived 
from  the  first  arrivals. 
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Figure  11.  Top  left:  waveform  just  before 
station  arrival  from  the  structure  with  the 
central  perturbation.  Top  right:  same  with  the 
offset  perturbation.  Right:  Hilbert  transformed 
cross-correlation  waveforms  showing  travel 
time  differences. 
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Measurements  of  the  travel  times  from  the  cross-correlated  waveforms  shown  in  Figure  1 1  are 
compared  with  the  predictions  from  ray  theory  and  from  using  the  kernel  above  for  the  three  cases 
in  Table  3.  The  offset  perturbation  does  cause  a  decrease  in  arrival  time  as  predicted  by  the  kernel, 
although  the  predicted  decrease  is  larger  than  measured.  For  the  central  and  large  perturbations. 
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ray  theory  is  considerably  more  accurate  than  the  kernel  integral.  This  is  probably  because  the 
Born  approximation  used  in  equation  4  is  valid  only  for  small  perturbations,  and  the  velocity 
change  here  is  large  enough  that  it  may  exceed  the  range  of  validity.  Note  that  for  ray  theory  in 
Table  3  we  have  used  the  exact  relation  for  travel  time:  AT  =  r/c-r/c0  instead  of  the  approximate 
equation  18,  which  gives  21.6  seconds,  a  result  more  similar  to  the  kernel  for  the  large  case. 


Table  3.  Measured  and  predicted  travel  time  variations  for  low  velocity  perturbations  (seconds) 


Perturbation 

Measured 

Ray  Theory 

Kernel 

Central 

6.63 

6.56 

2.72 

Offset 

-0.68 

0 

-1.25 

Large 

27.0 

26.26 

20.26 

4.3.3  Limitations  of  sensitivity  kernels 

We  investigate  the  extent  to  which  the  phase  delay  is  linearly  related  to  perturbation  in  phase 
velocity  through  the  sensitivity  kernels  calculated  using  the  membrane  model.  We  ran  three  test 
cases  with  perturbations  at  the  locations  shown  in  Figure  10.  In  each  case,  we  computed  the 
propagation  for  velocity  perturbations  of  -1,  -2,  -5,  -10,  and  -15%  in  a  background  phase  velocity 
of  3200  m/s.  (i.e.  low  velocity  regions  of  3168,  3136,  3040,  2880,  and  2720  m/s).  The  location  of 
the  perturbation  in  the  three  cases  correspond  to  1)  a  small  region  along  the  direct  path  between 
the  source  and  receiver;  2)  a  small  region  in  the  “first  Fresnel  zone”  predicted  by  the  kernel  to 
produce  a  phase  advance  from  the  direct  path;  and  3)  a  region  completely  crossing  the  kernel 
through  the  1st  Fresnel  zone. 

Waveforms  for  the  small  central  perturbation  are  shown  in  Figure  12.  As  expected,  the  phase  delay 
increases  with  the  increasingly  negative  velocity  perturbation.  Waveforms  for  the  offset 
perturbation  are  in  the  following  Figure  13.  The  phase  differences  in  this  second  case  are  subtle, 
but  the  modeling  results  in  a  waveform  that  has  the  predicted  phase  advance  for  a  negative  velocity 
perturbation.  Figure  14  and  Figure  15  show  snapshots  of  the  wave-field  at  the  time  of  arrival  at  the 
station  for  the  membrane  model  with  the  15%  negative  velocity  perturbation  centered  and  offset. 
In  both  cases  a  diffracted  waveform  is  apparent.  It  is  the  interference  of  the  diffracted  wave  and 
the  direct  wave  that  results  in  the  phase  advance  and  delay. 
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Figure  12.  Waveforms  for  the  1,2,5,10  and  15%  negative  perturbation  located  in  a  small  central 
area,  as  shown  in  (earlier)  Figure.  Top  panel  shows  the  calculated  waveforms;  the  bottom  panel 
shows  the  waveforms  correlated  with  the  unperturbed  signal.  A  consistent  phase  delay  is 
introduced  for  each  of  the  perturbation  cases,  along  with  a  signal  distortion. 


1  5 


Figure  13.  Waveforms  for  the  1,  2,  5, 10  and  15%  negative  perturbation  located  in  an  area  offset 
from  the  direct  ray  path.  Top  panel  shows  the  calculated  waveforms;  the  bottom  panel  shows  the 
waveforms  correlated  with  the  unperturbed  signal.  This  time,  a  phase  advance  appears  for  the 
perturbation  cases,  1,  2,  and  5%,  along  with  an  increasing  signal  distortion  caused  by  diffraction. 
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Figure  14.  Snapshots  of  wave  field  for  small  central  area  with  negative  velocity 
perturbation  of  15%.  The  source  is  located  at  x=100  km,  y=1200  km.  The  snapshot  is 
at  the  arrival  time  for  a  receiver  at  x=1700  km,  y=1200  km 
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Figure  15.  Snapshots  of  wavetleld  for  small  offset  area  with  negative  velocity 
perturbation  of  15%.  Snapshot  parameters  as  in  Figure  4. 

The  phase  delays  for  the  various  cases  are  presented  in  the  following  Table  4.  For  the  central 
perturbations,  ray  theory  and  the  model  agree  for  small  perturbations,  with  the  model  calculation 
showing  increasingly  faster  at  the  larger  perturbations.  This  agrees  with  the  idea  that  the 
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narrowband  membrane  waves  are  seeing  an  average  velocity  over  a  broad  region  outside  that  of 
the  perturbation. 

For  the  offset  perturbation,  the  kernel  predicted  phase  advance  for  a  negative  velocity  perturbation 
is  roughly  of  the  correct  size  for  the  smallest  perturbations.  However  for  larger  perturbations,  the 
predicted  phase  advance  does  not  agree  well  with  the  measured  phase  arrival  times. 


Table  4.  Phase  delays  for  kernel  and  direct  ray  test  cases. 


Perturbation 

1% 

2% 

5% 

10% 

15% 

Area 

Small  Central 
Model 

0.320 

0.625 

1.500 

3.024 

5.122 

Small  Central 
Kernel 

0.1898 

0.380 

0.949 

1.898 

2.847 

Small  Central; 
Direct  ray 

0.316 

0.638 

1.645 

3.472 

5.515 

Offset 

Model 

-0.094 

-0.132 

-0.060 

0.173 

0.324 

Offset 

Kernel 

-0.0798 

-0.160 

-0.399 

-0.798 

-1.197 

Large  Central 
Model 

1.318 

2.665 

6.878 

14.4 

22.4 

Large  Central 
Kernel 

1.226 

2.445 

6.113 

12.26 

18.339 

Large  Central 
Ray 

1.2626 

2.551 

6.579 

13.889 

22.059 
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4.4  Test  of  SurfMembrane  on  KNET  Structure 


Stevens  et  al  (2008)  performed  a  3D  calculation  of  propagation  of  seismic  waves  from  Lop  Nor  to 
KNET,  and  used  this  for  comparison  with  a  Bom  approximation  of  surface  waves  passing  through 
the  same  structure.  We  now  do  the  same  comparison  using  membrane  surface  waves. 

We  used  the  3D  elastic  finite  difference  code  TRES3D  to  calculate  wave  propagation  in  an  earth 
model  corresponding  to  the  region  between  Lop  Nor  and  KNET.  The  numerical  grid  size  is  678 
(North)  X  858  (East)  X  202  (Depth)  covering  an  area  of  approximately  (36-48N)  and  (71.5-92.5E). 
The  grid  spacing  in  each  dimension  is  2km  and  time  step  is  0.12s. 

The  heterogeneous  grid  model  is  generated  from  the  1°  x  1°  global  model  of  Stevens  et  al.  (2005) 
and  the  homogeneous  model  is  obtained  at  the  LopNor  (41.55N,88.70E)  source  site.  No 
attenuation  is  included  in  the  numerical  calculations.  The  homogenous  model  is  shown  in  Figure 
16  and  a  slice  through  the  heterogeneous  model  is  shown  in  Figure  17. 
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Figure  16.  Homogeneous  layered  model.  Numbers  show  Vp,  Vs  and  density.  Trailing  zeros  indicate 
infinite  Q. 
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Figure  17.  Heterogeneous  model  vertical  profile  of  Vp  (left)  and  Vs  (right)  at  about  42N. 
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Figure  18  shows  the  10  second  phase  velocity  map  calculated  from  this  model.  Figure  19  shows 
the  vertical  ground  motion  calculated  at  two  times  from  the  finite  difference  calculation,  and  Figure 
20  shows  the  10  second  membrane  surface  wave  calculation  at  the  same  times. 
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Figure  18. 10  second  phase  velocity  map  calculated  from  KNET  earth  model. 
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Figure  19.  3D  Unite  difference  calculation  of  vertical  ground  motion  from  an  explosion  at  Lop  Nor. 


Figure  20.  Membrane  surface  wave  calculated  at  120  seconds  and  240  seconds  using  the  10  second 
phase  velocity  map. 
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Figure  21  and  Figure  22  show  the  amplification  at  10  and  20  seconds  of  the  maximum  amplitude 
at  each  point  in  the  complex  structure  relative  to  the  maximum  amplitude  in  a  uniform  structure. 
The  finite  difference  and  membrane  surface  wave  calculations  give  very  similar  results,  both 
showing  strong  diffraction  by  heterogeneous  structures.  They  are  not  identical,  however,  since  the 
membrane  surface  wave  calculation  is  limited  to  10  seconds  period  while  the  finite  difference 
calculation  is  broadband.  Because  the  variations  are  due  to  diffraction,  small  differences  in 
frequency  band  can  change  the  diffraction  patterns  substantially.  Our  main  interest  is  in  assessing 
whether  membrane  surface  waves  provide  a  practical  technique  for  calculating  amplitude 
variations  in  heterogeneous  structure,  so  we  look  at  that  question  in  more  detail. 


200  400  600  800  1000  1200  1400  1600  200  400  600  800  1000  1200  1400  1600 


Figure  21.  ln(amplification)  at  10  seconds.  Left:  FD,  Right:  membrane.  “Amplification”  is  the  ratio 
of  amplitude  in  the  heterogeneous  structure  to  the  homogeneous  structure. 


We  calculated  the  ratio  of  the  maximum  amplitude  of  the  waves  in  heterogeneous  to  homogeneous 
structure,  band-pass  filtering  close  to  10  seconds  period  (0.09  to  0.11  Hz)  for  all  points  in  the 
calculation.  An  example  of  the  filtered  waveforms  is  shown  in  Figure  23.  We  want  to  assess  how 
well  the  membrane  surface  wave  calculations  can  be  used  to  correct  amplitude  variations  from  a 
3D  structure.  That  is,  we  would  like  to  be  able  to  calculate  the  amplification  using  SurfMembrane 
to  generate  an  amplitude  correction  at  each  point  that  can  be  subtracted  from  the  data  (in  this  case 
the  3D  calculation)  to  produce  the  amplitudes  for  the  same  event  in  a  homogeneous  structure. 
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Figure  23.  Filtered  waveforms  at  a  single  point  for  the  heterogeneous  (top)  and  homogeneous 
(bottom)  structures,  calculated  with  membrane  surface  waves  (left)  and  3D  finite  difference  (right). 
Both  show  a  similar  amplitude  reduction  in  the  complex  structure. 

We  used  the  membrane  surface  wave  calculation  to  “correct”  the  heterogeneous  finite  difference 
calculation.  For  each  point  we  calculated  the  ratio  of  the  heterogeneous  to  homogeneous 
membrane  surface  wave  calculation,  and  divided  the  heterogeneous  finite  difference  ratio  by  this 
value.  If  the  correction  were  perfect  and  the  variation  generated  by  SurfMembrane  exactly  equaled 
the  variation  in  the  3D  finite  difference  calculation,  then  this  ratio  would  be  one  at  every  point. 
The  results  displayed  as  a  histogram  are  shown  in  Figure  24.  While  not  every  point  is  improved, 
the  amplitude  distribution  is  substantially  narrowed,  and  the  correction  appears  to  be  particularly 
effective  for  low  amplitudes,  which  are  largely  removed  in  the  corrected  data  set.  The  standard 
deviation  in  the  amplitude  ratios  is  reduced  from  0.22  to  0.15.  This  could  probably  be  improved 
by  varying  the  filter  parameters,  and  this  shows  that  if  the  earth  structure  is  known  well  enough, 
the  membrane  calculation  can  be  used  to  correct  surface  wave  amplitudes. 
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Figure  24.  Histograms  of  log  amplitude  distribution  in  the  finite  difference  calculation  (left),  and  in 
the  finite  difference  calculation  corrected  with  the  membrane  surface  wave  calculation  (see  text). 
The  spread  is  substantially  reduced,  and  low  amplitudes  in  particular  are  greatly  improved. 
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4.5  SurfMembrane  Simulation  of  US  Array  Data 

SurfMembrane  was  run  for  all  73  events  recorded  by  the  US  Array  discussed  earlier.  All  were  run 
with  the  Leidos  models,  and  events  in  the  western  US  were  also  run  with  the  CU  models.  An 
example  is  shown  in  Figure  25  for  the  2005/06/12  M  5.2  Anza  earthquake,  for  10  second  period 
waves.  Both  calculations  show  effects  of  diffraction  during  propagation.  The  CU  model  does  not 
include  the  ocean,  so  the  effects  of  reflection  and  scattering  from  the  Ocean/continent  boundary 
are  seen  only  in  the  Leidos  model. 
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Figure  25.  SurfMembrane  calculation  of  waveforms  from  the  Anza  earthquake  propagating 
through  the  CU  (left)  and  Leidos  (right)  models. 

Path  dependent  station  corrections  were  derived  for  all  events  described  in  section  3  that  were 
recorded  on  US  Array  stations.  The  purpose  of  the  station  corrections  is  to  remove  the  effects  of 
the  heterogeneous  structure  so  that  the  amplitudes  correspond  to  a  common  uniform  structure. 
Station  corrections  were  derived  using  two  different  sets  of  sources.  In  both  cases,  the  approach 
was  to  run  SurfMembrane  for  the  heterogeneous  structure  and  a  homogeneous  structure,  and  then 
subtract  the  logarithm  of  the  homogenous  amplitude  from  the  logarithm  of  the  heterogeneous 
amplitude.  This  gives  a  station  residual  for  the  effect  of  heterogeneity  that  can  be  used  as  a  station 
correction.  The  first  set  of  sources  corresponded  to  the  known  CMT  solutions  for  each  source.  The 
second  set  used  point  explosion  sources  at  0.5  km  depth.  Although  the  CMT  sources  should  more 
accurately  reproduce  the  observed  signals  from  each  event,  they  have  two  disadvantages:  1)  in  the 
general  case  when  we  want  to  apply  station  corrections,  in  most  cases  we  will  not  have  a  CMT 
solution;  and  2)  station  corrections  may  be  anomalously  large  or  small  near  nodes.  The  station 
corrections  for  all  paths  are  shown  in  Figure  26. 
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10  Second  Station  Corrections 
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Figure  26.  Ms  station  corrections  derived  from  Leidos  and  WUS  structures  for  all  events  analyzed 
that  were  recorded  on  US  Array  stations.  Left:  10  seconds;  right:  20  seconds.  Top:  CMT  sources. 
Bottom:  Explosion  sources.  Variation  is  larger  for  CMT  sources  because  differences  are  amplified 
near  nodes. 


We  expect  the  explosion  sources  to  generate  the  same  focusing/defocusing  and  diffraction  effects 
as  the  CMT  sources  and  to  be  more  robust  against  nodal  anomalies,  so  in  general  we  expect  them 
to  be  the  better  choice  for  application  of  station  corrections.  The  mean  and  standard  deviation  of 
the  station  corrections  as  a  function  of  range  is  shown  in  Figure  27.  As  can  be  seen  in  both  Figure 
26  and  Figure  27  there  is  more  variability  at  10  seconds  than  at  20.  There  is  also  more  variability 
in  the  station  corrections  derived  from  the  Leidos  structures  than  in  the  station  corrections  derived 
from  the  WUS  structures,  most  likely  due  to  the  Leidos  structures  having  sharper  boundaries. 
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Figure  27.  Mean  and  standard  deviation  of  station  corrections  as  a  function  of  range  for  Leidos 
structures  (top)  and  WUS  structures  (bottom).  Left:  10s;  right:  20s. 

Figure  28  and  Figure  29  show  comparisons  of  observed  and  calculated  amplitudes  for  6  events  at 
20  and  10  seconds  period,  respectively.  Calculations  are  shown  for  both  the  WUS  and  Leidos 
structures.  Although  the  calculated  and  observed  amplitudes  agree  fairly  well,  it  is  clear  that  almost 
all  of  the  amplitude  variation  is  due  to  radiation  pattern  rather  than  structural  effects.  Furthermore, 
the  comparison  is  very  sensitive  to  accurate  CMT  solutions.  For  event  37  in  particular,  the 
difference  between  observed  and  calculated  amplitudes  is  very  large  at  some  azimuths  due  to  errors 
in  the  predicted  radiation  pattern.  The  observed  and  calculated  data  have  very  similar  shapes  for 
that  event,  but  they  are  shifted  by  about  30  degrees.  Both  the  calculated  and  observed  variability 
due  to  structure  are  small  at  20  seconds  (remember  that  all  paths  are  less  than  1000  km). 
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Figure  28.  Calculated  and  observed  20  second  surface  wave  amplitudes  for  6  events  recorded  on  the 
US  Array.  Black:  observations;  Blue:  calculations  with  WUS  %°  structures;  Magenta:  calculations 
with  Leidos  1°  structures . 
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Figure  29.  Calculated  and  observed  10  second  surface  wave  amplitudes  for  6  events  recorded  on  the 
US  Array.  Black:  observations;  Blue:  calculations  with  WUS  %°  structures;  Magenta:  calculations 


with  Leidos  1°  structures. 
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There  is  considerably  more  structural  variability  at  10  seconds,  and  the  amount  of  variability  in 
the  observed  and  calculated  data  is  comparable,  however  there  is  not  good  agreement  between 
calculated  and  observed  amplitude  variations  at  individual  stations. 

Figure  30  and  Figure  31  show  histograms  of  surface  wave  amplitudes  at  20  and  10  seconds, 
respectively,  for  one  event.  Two  types  of  corrections  are  used,  derived  from  SurfMembrane 
calculations  using  CMT  and  explosion  sources  as  described  above.  We  cannot  expect  a  large 
reduction  in  amplitude  variability,  because  most  of  the  amplitude  variability  is  due  to  radiation 
pattern.  With  the  CMT-derived  corrections  the  variability  is  actually  increased.  With  the 
explosion-derived  corrections  we  do  get  a  small  improvement.  At  20  seconds  the  standard 
deviation  is  reduced  from  0.280  to  0.266,  and  at  10  seconds  from  0.187  to  0.181.  The  10  second 
variability  is  smaller  than  the  20  second  variability  because  structural  effects  wash  out  some  of  the 
highs  and  (particularly)  lows  in  the  radiation  pattern  in  both  the  calculated  and  observed  data. 


Figure  30.  Histogram  of  20  second  surface  wave  amplitudes  for  one  event.  Left:  Uncorrected; 
Middle:  Corrected  with  CMT  derived  correction;  Right:  Corrected  with  explosion  derived 
correction. 


Figure  31.  Histogram  of  10  second  surface  wave  amplitudes  for  one  event.  Left:  Uncorrected; 
Middle:  Corrected  with  CMT  derived  correction;  Right:  Corrected  with  explosion  derived 
correction. 

For  some  events  and  some  stations  there  is  additional  variability  in  the  data  that  is  not  easily 
explained  by  either  radiation  pattern  or  structure.  Figure  32  shows  the  measured  amplitudes 
(shown  as  magnitudes  calculated  from  the  amplitudes)  for  one  event  and  the  amplitudes  calculated 
with  the  two  sets  of  structures.  The  calculations  reproduce  the  observed  radiation  pattern  quite 
well,  especially  at  20  seconds,  however  there  is  still  a  substantial  amount  of  unexplained  variability 
in  the  data  including  high  observed  amplitudes  right  at  the  node. 
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We  took  a  closer  look  at  the  station  with  the  highest  amplitude,  which  is  located  very  close  to  the 
node.  It  would  take  an  extreme  structural  variation  to  cause  this  type  of  effect.  That  station  is 
H17A,  which  is  located  just  west  of  Yellowstone  Lake  and  event  21  is  just  west  of  Reno,  NV.  The 
path  does  have  some  significant  complexity,  but  the  earth  model  we  are  using  does  not  have 
enough  variability  to  cause  an  effect  of  that  magnitude.  However,  Lou  and  van  der  Lee  (2014) 
show  a  very  large  S-wave  travel  time  anomaly  at  exactly  the  same  location,  so  there  may  be  a  large 
structural  anomaly  that  is  not  captured  by  our  current  models. 


Event  21  10  sec 


Azimuth  (deg) 


Figure  32.  Predicted  and  observed  station  magnitudes  for  event  21.  Black  crosses  are  observed 
magnitudes,  blue  x  are  WUS  calculations,  magenta  x  are  Leidos  calculations.  Calculations  used  the 
CMT  solution  for  each  source.  Left:  10s,  right:  20s. 

Figure  33  shows  amplitude  vs.  distance  for  stations  close  to  the  direction  of  H17A  and  although 
there  is  some  increase  in  amplitude  for  the  two  next  most  distant  stations,  this  station  is  clearly  an 
outlier.  Figure  33  also  shows  waveforms  from  H17A  and  I16A,  which  is  in  the  same  direction,  but 
100  km  closer  to  the  source  and  has  only  20%  of  the  amplitude.  Both  waveforms  look  normal  and 
show  no  evidence  of  data  problems  such  as  a  noise  spike.  It  is  possible  there  is  a  calibration  error 
at  this  station,  however  a  comparison  of  several  other  events  also  recorded  at  this  station  show 
considerable  variation  but  not  a  persistent  high  anomaly. 
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Figure  33.  Top:  amplitudes  at  7  stations  close  to  50  degree  azimuth.  Bottom:  waveforms  at  station 
I16A  (left)  at  850  km  and  H17A  (right)  at  942  km  band-pass  filtered  at  10  seconds.  The  amplitude 
increases  by  a  factor  of  5  between  the  two  stations. 
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4.6  Inversion  of  Eurasian  Attenuation  Data  for  Q  Structure 


Inversion  of  the  Eurasian  attenuation  data  set  follows  the  same  procedure  used  by  Stevens  et  al 
(2008),  which  is  summarized  here.  Inversion  of  attenuation  data  for  Q  structure  is  accomplished 
by  solving  equation  19  using  the  LSQR  algorithm  (Nolet,  1987). 


where: 


1.  The  data  Ad  are  attenuation  residuals.  Attenuation  estimates  are  derived  from  an 
existing  Q  model,  and  the  differences  between  those  and  the  observations  are  the  data 
used  in  the  inversion. 


2.  The  matrix  A  is  derived  from  partial  derivatives  of  the  attenuation  coefficients  with 
respect  to  model  Q  in  each  layer  for  the  path-averaged  inversion.  The  velocities  and 
densities  used  to  calculate  partial  derivatives  correspond  to  the  final  models  from 
Stevens  et  al  (2005). 

3.  The  starting  model  x  and  constraining  model  xc  are  the  same,  and  are  derived  from 
PREM  for  depths  greater  than  100  km,  and  Swanger’s  Law  (Q=100  p,  with  p  in  km/s) 
for  shallower  depths.  The  values  derived  from  PREM  are  Q=18  P  for  depths  between 
100  km  and  220  km,  Q=30  p  at  greater  depths.  There  are  discontinuities  at  100  km 
and  220  km  and  a  smoothness  criterion  is  applied  to  layers  above  and  below  100  km. 
The  inversion  is  performed  for  layers  shallower  than  220  km.  Q  is  fixed  to  1 8  p  at 
220  km  depth  and  to  30  P  below  this  depth. 

4.  H  is  a  difference  operator  that  applies  to  vertically  neighboring  layers  and  has  the 
effect  of  constraining  the  vertical  smoothness  of  Q  profile.  H  applies  to  layers  in  the 
crust  and  upper  mantle,  but  has  explicit  discontinuities  as  described  above,  s  is  the 
weighting  of  the  smoothness  constraint  and  can  be  a  diagonal  matrix  (for  variably 
weighted  smoothing)  or  a  scalar  (constant  smoothing).  I  is  the  identity  matrix  and  X 
weights  the  damping  which  constrains  the  norm  of  the  difference  between  final  Q 
structures  and  constraining  model  xc.  X  can  be  a  scalar  for  constant  damping,  or  a 
diagonal  matrix  for  variable  damping. 


5.  The  model  vector  Ax  consists  of  p/Q  for  each  layer  that  is  free  to  change  in  each 
structure,  and  optionally  can  include  the  change  in  moment  of  each  event.  That  is,  for 
attenuation  residuals  that  were  derived  using  spectra  with  a  fixed  model  moment,  the 
moment  for  each  event  can  be  corrected  as  part  of  the  inversion. 

6.  8  is  the  vector  of  residuals  that  remain  after  inversion  (the  inversion  minimizes  le|). 

Unlike  inversion  for  shear  velocity,  inversion  for  Q  is  linear,  so  only  a  single  iteration  is  necessary, 
although  multiple  inversions  are  done  with  different  damping  and  smoothing  parameters  to 
generate  realistic  models. 

Surface  wave  propagation  in  a  heterogeneous,  anelastic  structure  takes  the  following  form, 
separating  source,  path  and  receiver  (notation  follows  Harkrider  et  al,  1994): 
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u  A(0,r,(p)  = 


2  A, 


\[ae  sin(r/^) 


ncoc. 


T-\/c2AR2  exP  (V 4  -  ^  /  cp  -  ypr)]  Fs  (0,  q>,  h ) 
1 


(20) 


where  0  is  angular  frequency,  r  is  source  to  receiver  distance,  h  is  source  depth,  ae  is  the  radius  of 
the  earth,  <p  is  azimuth,  Ar  is  the  Rayleigh  wave  amplitude  function,  c  is  phase  velocity,  y  is  the 
attenuation  coefficient,  and  the  subscripts  1,  2,  and  p  refer  to  parameters  derived  from  the  source 
region  structure,  parameters  derived  from  the  receiver  region  structure,  and  parameters  which  are 
defined  by  path  averages,  respectively.  All  source  properties  are  contained  in  the  function  Fs. 

Collecting  the  distance  independent  terms,  equation  20  can  be  written  in  the  form: 


uz(co,  r)\  =  M0 


S(0)A(0)exp[-yp(0)r] 
■s]ae  sin(>  /  ae) 


(21) 


where  Mo  is  the  source  moment  and  S(co)  the  source  function  (important  for  large  events),  A(co)  is 
a  frequency  dependent  function  that  depends  on  source  and  receiver  structure  and  focal 
mechanism.  We  assume  that  A(co)  can  be  predicted  well  enough  from  the  background  earth 
structure  and  the  source  mechanism  which  is  usually  defined  by  a  CMT  moment  tensor.  While  the 
inversion  program  has  the  capability  to  allow  an  amplitude  scale  factor  for  each  path,  which  would 
allow  variations  in  explosion  amplitude  due  to  tectonic  release,  and  variations  in  earthquake 
amplitude  due  to  errors  in  source  mechanism,  in  the  inversions  that  follow  we  only  allowed  the 
moment  to  vary,  which  is  a  constant  factor  for  all  paths  for  a  single  event.  In  equation  21,  Mo  and 
yp  are  derived  from  a  starting  source  mechanism  and  background  earth  model,  and  allowed  to  vary 
in  the  inversion  while  the  other  factors  are  held  fixed.  S(co)  is  derived  for  a  triangular  function  with 
rise  time  (half-duration)  T.  Since  this  is  approximate,  points  where /7>0.5,  which  corresponds  to 
an  amplitude  reduction  of  0.4,  are  zero  weighted.  The  observed  data  can  then  be  written: 


I  or  J  t/ro  S(o))A(o))exp[-y  (0)r] 

K  (0,  r)  =  M  °0 -  ===== - 

y]aesm(r/ae) 


(22) 


And  the  log  ratio  of  observed  to  predicted  spectra  has  the  form: 


In 


u°(co,  r ) 


u(co,r ) 


-  In  + 1 (, o))r  -  yp°  (0)7-]  =  A  In  M0  +  rAy 


M 


(23) 


Ay  can  be  further  expanded  into  a  sum  over  Q  structure  in  each  structure  traversed  by  the  ray  along 
the  source  to  receiver  point  multiplied  by  the  fraction  of  ray  over  each  structure.  So  we  can  rewrite 
equation  23  as: 


u°(a),r) 
uz(a>,  r) 


1=1 


rik  7=1 


(24) 


The  subscript  k  refers  to  an  event,  i  refers  to  a  single  path  for  event  k,  1  refers  to  each  model  type 
traversed  by  the  ray  and  Arljk  / rlk  is  the  fraction  of  the  path  that  the  ray  spends  in  each  model.  The 

subscript  j  refers  to  layer  number  in  each  model.  L  is  the  total  number  of  model  types  traversed 
and  J  is  the  total  number  of  layers  allowed  to  change  in  each  model.  So  the  data  in  the  inversion  is 
the  left  hand  side  of  equation  24,  the  spectral  ratio  divided  by  the  distance  for  multiple  frequencies, 


Approved  for  public  release;  distribution  is  unlimited. 

32 


and  the  inversion  is  performed  for  the  quantities  A  In  M * ,  the  change  in  moment  of  each  event,  and 
A  (j3/Q),.,  the  ratio  of  shear  velocity  to  Q  in  each  layer  of  each  model.  The  function  Gij  gives  the 

change  in  y  with  respect  to  change  in  p/Q  and  can  be  written  assuming  no  bulk  attenuation  in  terms 
of  partial  derivatives  of  phase  velocity  with  respect  to  material  velocities  as: 

Gi=Al +iAA_  (25) 

dptj  3  au  dav 


The  inversions  described  below  all  invert  data  of  the  form  described  above.  We  also  compare  with 
interstation  attenuation  estimates.  Using  equation  23  again  at  two  or  more  stations  at  distances  rn, 
we  get  a  set  of  equations  of  the  form: 


In 


uz(0J,rn) 


=  AlnM0  +  rnAy 


(26) 


The  equations  define  a  line  with  slope  Ay  and  intercept  AlnMo  so  this  can  be  used  as  a  check  on 
the  inversion  by  finding  the  average  change  in  attenuation  over  all  or  a  subset  of  paths  for  each 
event,  and  the  change  in  moment. 

4.6.1  Data  used  in  the  Q  Inversion 

As  discussed  in  section  3,  we  retrieved  59,000  waveforms  from  1850  Eurasian  and  African  events. 
We  also  retrieved  the  locations,  depths  and  origin  times  for  all  of  these  events  from  the 
International  Seismological  Center  (ISC),  and  replaced  the  CMT  origins  and  depths  with  the  better 
ISC  origins.  We  then  derived  phase-matched  filters  from  the  CMT  moment  tensors  of  each  event, 
applied  these  to  the  data,  and  assembled  a  data  set  of  smoothed,  instrument  corrected  spectra. 
These  measured  spectra  were  then  divided  by  the  theoretical  calculated  spectral  amplitudes  and 
saved  at  a  set  of  frequencies  corresponding  to  periods  between  8  and  50  seconds. 
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Figure  34.  Left:  Map  showing  location  of  all  rays  used  in  the  inversion.  Right:  histogram  of  path 
lengths  used  in  the  inversion. 


These  formed  the  data  on  the  left  hand  side  of  equation  26.  These  data  were  then  combined  with 
the  data  set  from  the  2008  study  and  put  through  an  aggressive  data  quality  procedure  to  prune  out 
bad  or  noisy  data,  and  deep  events  or  data  with  bad  depths  or  focal  mechanisms.  This  resulted  in 
a  high  quality  data  set  of  23,148  spectra  from  998  events.  There  were  an  average  of  14.8  spectral 
data  points  for  each  waveform,  leading  to  a  total  data  set  of  343,322  data  points.  Figure  34  shows 
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a  map  of  all  of  the  ray  paths  used  in  the  inversion,  and  the  distribution  of  path  lengths  in  the  data 
set. 

4.6.2  Q  Inversion  Results 

The  result  of  the  inversion  is  adjustments  to  Q  structure  in  each  model  layer  and  adjustments  to 
the  moment  of  each  event.  Figure  35  shows  the  adjustment  in  moments  (equation  26)  that  were 
output  from  the  inversion.  Figure  35  also  shows  the  mean  and  mean  ±  standard  deviation  of  the 
raw  data,  the  data  after  weighting  (derived  from  quality  control  procedure)  was  applied,  and  after 
moment  adjustments  were  applied. 
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Figure  35.  Left:  Adjustments  to  moments  in  the  CMT  data  set.  Right:  average  gamma  data  over  all 
paths,  together  with  average  gamma  after  moment  adjustments  and  weighting. 


Figure  36  shows  a  comparison  of  inversion  results  with  the  data  and  with  the  starting  model.  As 
we  also  found  in  the  2008  study,  the  data  and  inversion  results  show  considerably  more  attenuation 
than  the  starting  model,  and  much  more  variability.  The  standard  deviations  are  larger  at  the  closer 
ranges  because  variations  due  to  structure  are  larger  relative  to  variations  due  to  attenuation  at 
short  distances. 
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Figure  36.  Gamma  data  and  inversion  results  compared  with  starting  model  over  several  distance 
ranges.  Upper  left:  all  data;  upper  right:  300-2000  km;  lower  left:  1000-5000  km;  lower  right:  5000- 
10000  km.  Solid  lines  are  mean;  dashed  lines  mean  ±  standard  deviation.  Numbers  in  the  titles  are 
distance  ranges. 


Figure  37  shows  the  results  of  the  inversion  plotted  on  a  map  of  Eurasia  and  Africa.  There  is  a 
band  of  high  surface  wave  attenuation  that  stretches  across  the  Middle  East  from  the  Mediterranean 
Sea  to  India.  This  is  also  consistent  with  the  2008  results,  although  there  is  much  more  data  in  this 
inversion  than  the  2008  inversion. 
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Figure  37.  Left:  Attenuation  maps  at  four  frequencies.  From  top  to  bottom:  0.02  Hz  (50  seconds 
period),  0.05  Hz  (20  seconds),  0.067  Hz  (15  seconds),  0.1  Hz  (10  seconds).  Right:  histograms  of 
gamma  values  for  each  frequency. 
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4.7  Application  of  SurfMembrane  to  Eurasian  events 


Surface  wave  membrane  calculations 
require  considerably  more  computation 
time  for  the  Eurasian  events  because 
paths  are  much  longer,  averaging  about 
5,000  km.  Calculations  take 
approximately  one  hour  to  complete, 
still  much  less  than  a  corresponding  3D 
finite  difference  calculation.  For  these 
calculations,  we  derived  phase  velocity 
maps  centered  on  each  event  (Figure 
38),  and  we  ran  four  calculations  for 
each  event:  a  calculation  using  the 
CMT  solution  with  the  heterogeneous 
structure  and  a  calculation  using  the 
CMT  solution  with  a  homogeneous 
structure,  and  explosion  calculations  in 
the  heterogeneous  and  homogenous 
structures. 
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Figure  38.  Example  phase  velocity  map  for  event  20300. 


Membrane  calculations  in  heterogeneous  structures  generated  surface  waves  with  much  more 
variability  than  uniform  structures.  Figure  39  shows  the  surface  waves  for  two  calculations  at  10 
and  20  second  periods  after  650  time  steps  at  which  time  they  have  propagated  to  about  4000  km 
from  the  source.  All  show  strong  signs  of  diffraction  and  the  surface  waves  seem  to  break  up  where 
they  cross  strong  boundaries  such  as  at  boundaries  with  oceanic  structures,  particularly  at  10 
seconds.  Figure  40  shows  the  peak  amplitudes  at  each  point  from  the  homogeneous  and 
heterogeneous  calculations  for  one  of  these  events. 
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Figure  39.  Surface  waves  from  two  events  at  10  seconds  (top)  and  20  seconds  (bottom)  showing 
strong  diffraction,  especially  for  10  seconds  waves  crossing  oceanic  boundaries. 
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Figure  40.  Log  surface  wave  peak  amplitudes  at  each  point  for  event  20300  at  10  seconds  period 
(left)  and  20  seconds  (right). 


Figure  41  shows  the  ratio  of  the  peak  amplitudes  between  the  heterogeneous  and  homogeneous 
calculations  at  each  point.  These  “dandelion  plots”  correspond  to  the  predicted  surface  wave 
amplitude  correction  for  each  point,  and  they  can  be  quite  large.  Unfortunately,  they  are  also  very 
sensitive  to  earth  structure  along  the  path  and  small  changes  in  earth  structure  can  cause  these 
patterns  to  shift,  substantially  changing  the  amplitude  correction  at  each  point. 
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Figure  41.  Log  ratio  of  surface  wave  amplitudes  in  heterogeneous  to  homogeneous  structure  for  the 
two  events  shown  in  Figure  39.  Top:  10  seconds.  Bottom:  20  seconds. 


4.7.1  Ms  Estimation  in  Heterogeneous  Structures 

We  want  to  use  surface  wave  membrane  calculations  to  correct  for  structural  variations  in 
amplitude.  In  this  section  we  perform  an  uncertainty  analysis  to  determine  under  what  conditions 
such  corrections  are  likely  to  improve  Ms  measurements.  The  membrane  calculations  rely  on  phase 
velocity  maps  acquired  through  tomographic  inversion.  The  phase  velocity  maps  are  then  used  to 
calculate  the  effects  of  material  properties  on  surface  wave  propagation.  While  the  calculations 
can  model  diffraction  and  scattering  in  heterogeneous  structures,  the  results  are  only  as  accurate 
as  the  underlying  Earth  model.  Any  errors  in  the  phase  velocity  map  used  to  perform  the 
calculation  are  propagated  into  the  subsequent  surface  wave  magnitude  estimate.  In  the  analysis 
above,  calculation  of  “dandelion”  plots  that  estimate  the  Ms  corrections  from  the  surface 
membrane  calculations  have  assumed  that  the  underlying  models  are  approximately  valid.  Here 
we  examine  the  effects  of  phase  velocity  model  uncertainties  on  the  resulting  surface  waves. 
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Simulation  Approach  &  Uncertainty  Estimation 


Figure  42.  Source  locations  for  the  88  surface  membrane  calculations  distributed  throughout 
Eurasia  and  parts  of  Africa.  The  numbers  near  the  red  squares  are  an  identifying  label. 


For  the  following  analysis,  we  use  SurfMembrane  to  calculate  surface  wave  magnitudes  at  20- 
second  period.  Although  surface  wave  measurements  are  performed  over  a  narrow  band  of 
frequencies,  filtering  still  contains  components  from  nearby  frequencies,  so  we  must  account  for 
the  subtle  differences  in  phase  velocity  by  generating  results  by  superposition  of  at  least  a  few 
phase  velocities.  To  approximately  account  for  the  frequency  dependence  of  phase  velocity,  we 
combine  the  results  of  three  separate  surface  membrane  calculations  at  phase  velocities  of  19,  20, 
and  21  seconds.  Eighty-eight  source  locations  distributed  throughout  Eurasia  and  parts  of  Africa 
were  tested  using  an  explosion  source  (Figure  42).  Each  calculation  is  centered  on  the  source 
location,  and  spans  ±3750  kilometers  distance  along  the  North-South  and  East-West  directions 
using  2  kilometer  grid  size  and  0.3  seconds  per  time  increment.  For  the  calculations  at  each  source 
location,  the  nominal  phase  velocity  map  is  obtained  by  projecting  the  global  Leidos  1 -degree 
model  to  a  two-dimensional  surface  centered  on  the  source  location  (Figure  43). 
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Figure  43.  Left  panel:  Phase  velocity  map  on  a  1-degree  grid  of  longitude  and  latitude.  Right  panel: 
Projection  of  the  map  in  the  left  panel  to  a  Northing-Easting  (x-y)  coordinate  system  using  the  local 
distance/bearing. 
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To  study  the  effects  of  uncertainties  in  the  phase  velocity  maps  on  the  calculated  surface  waves, 
perturbations  that  are  spatially  correlated  over  -100  kilometers  and  with  root-mean-squared 
fractional  amplitude  of  2.6%  are  added  to  the  nominal  phase  velocity  maps.  The  spatial  scale  and 
amplitude  of  the  perturbations  were  estimated  by  comparing  phase  velocity  maps  for  the  Western 
United  States  from  the  University  of  Colorado  (Shen  et  al.,  2013)  to  the  1-degree  maps  developed 
by  Leidos  (Stevens  et  al.,  2005)  in  the  same  region  (Figure  44). 


19  s  20  s  21s 


Fractional  Difference  in  Phase  Velocity 

Figure  44.  Top  panels:  Fractional  difference  in  the  WUSA  and  Leidos  phase  velocity  maps  at  19, 20, 
and  21  s  periods.  Bottom  panels:  Distribution  of  differences,  which  are  approximately  Gaussian 
with  standard  deviation  of  about  2.6%. 
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For  each  source  location,  the  nominal  Leidos  phase  velocity  model  is  assumed  to  be  the  “correct” 
phase  velocity  map.  An  ensemble  of  twenty-five  instances  of  perturbation  maps  is  applied  to  the 
nominal  model,  and  the  surface  membrane  calculation  is  run  for  each  instance.  An  example  of  an 
instance  of  the  perturbation  map  is  shown  in  Figure  45.  The  same  perturbation  map  is  applied  to 
each  phase  velocity  because  the  comparison  of  the  maps  in  Figure  44  shows  that  the  differences 
are  consistent  across  nearby  surface  wave  periods. 


Avg.  Fractional  Diff.  17-23  sec  (WUSA  vs.  Leidos)  Spatially  Correlated  Perturbation  Example 


Figure  45.  Left  panel:  Average  fractional  difference  in  the  WUSA  and  Leidos  phase  velocity  maps 
in  the  17-23  s  band.  Right  panel:  Example  of  phase  velocity  map  fractional  perturbation  that  is 
applied  to  the  Leidos  Earth  model. 

Each  perturbed  map,  when  added  to  the  nominal  map,  is  an  instance  that  theoretically  could  have 
been  obtained  given  other  data  used  in  the  inversion  (Figure  46).  When  the  underlying  phase 
velocity  map  changes,  the  resulting  surface  wave  amplitudes  also  change  (Figure  47).  Macroscopic 
features  of  the  surface  wave  amplitudes,  such  as  the  number  and  approximate  direction  of 
amplification  zones,  are  usually  preserved.  However,  the  precise  directions  where  amplification 
reaches  a  local  maximum  or  minimum  can  shift  by  a  few  degrees.  At  long  ranges  even  a  small 
shift  in  direction  will  move  an  amplification  zone  hundreds  of  kilometers.  In  other  words, 
uncertainty  in  the  phase  velocity  of  even  a  few  percent  (due  to  limited  data,  systematic  errors,  or 
for  any  other  potential  reason)  can  move  a  particular  location  from  a  region  of  amplification  to  a 
region  of  de-amplification.  This  shows  the  importance  of  minimizing  errors  and  uncertainty  in  the 
generation  of  phase  velocity  maps. 


Instance  1 


Instance  2 


Instance  3 


Figure  46.  Examples  of  three  instances  of  perturbed  phase  velocity  maps  at  20-second  period.  The 
area  shown  is  a  7500  km  x  7500  km  region  near  the  Tarim  Basin ,  centered  on  source  #54.  Dark  blue 
color  is  2800  m/s  phase  velocity ,  and  dark  red  color  is  4200  m/s  phase  velocity. 
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Figure  47.  Examples  of  surface  wave  amplitudes  (dB)  produced  by  separate  perturbation  instances 
applied  to  the  same  nominal  phase  velocity  map.  The  area  shown  is  a  7500  km  x  7500  km.  The  color 
scale  spans  2.2  dBfrom  high  (red)  to  low  (blue)  amplification  regions.  Note,  these  surface  wave 
amplitudes  do  not  correspond  to  the  perturbation  instances  shown  in  Figure  6. 

For  each  of  the  88  source  locations,  the  root-mean-square  deviation  in  the  Ms  measurement  over 
the  ensemble  of  25  perturbation  instances  is  used  to  estimate  the  uncertainty  in  Ms  due  to 
uncertainty  in  the  underlying  phase  velocity  map.  Figure  8  shows  a  representative  example  of  an 
Ms  uncertainty  map. 
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Figure  48.  Ms  uncertainty  map  generated  by  calculating  the  standard  deviation  of  Ms  over  an 
ensemble  of  25  perturbation  instances.  This  is  from  the  same  source  as  shown  in  the  example  in 
Figure  47. 


Surface  Wave  Magnitude  &  Uncertainty  Analysis 

The  “correction”  to  the  Ms  estimate  due  to  heterogeneous  Earth  structure  is  obtained  by  subtracting 
the  surface  wave  magnitude  assuming  a  homogeneous  structure  from  the  magnitude  generated  by 
heterogeneous  structure.  Figure  49,  Figure  50  and  Figure  51  show  the  Ms  corrections  and  projected 
phase  velocity  maps  for  the  88  source  locations  throughout  Eurasia  and  parts  of  Africa  that  were 
shown  in  Figure  42.  Nearby  numbers  correspond  to  source  locations  in  proximity  to  each  other 
(see  Figure  42). 
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Figure  49.  Ms  corrections  (bottom)  and  projected  nominal  20-second  phase  velocity  maps  (top)  for 
source  locations  1-30. 
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Figure  50.  Ms  corrections  (bottom)  and  projected  nominal  20-second  phase  velocity  maps  (top)  for 
source  locations  31-60. 
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Figure  51.  Ms  corrections  (bottom)  and  projected  nominal  20-second  phase  velocity  maps  (top)  for 
source  locations  61-88. 
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The  highly  nonlinear  nature  of  the  resulting  surface  wave  magnitudes  makes  comparison  of  the 
entire  population  of  88  source  locations  difficult.  One  approach  taken  here  is  to  characterize  each 
source  location  by  the  angle-averaged  root-mean- square  (rms)  amplitude  of  the  Ms  corrections 
(restricted  to  ranges  less  than  2000  km,  as  explained  below).  For  a  source  embedded  in  a 
completely  uniform  structure  with  a  constant  phase  velocity,  the  Ms  correction  at  all  angles  and 
ranges  would  be  zero.  Thus,  the  rms  amplitude  of  the  Ms  correction  is  likely  associated  with 
gradients  in  the  phase  velocity  map. 

Figure  52  shows  the  gradient  of  the  20-second  phase  velocity  map  spanning  the  88  source 
locations.  The  rms  amplitude  of  the  Ms  correction  associated  with  each  source  location  is 
proportional  to  the  overlaid  scatter  point.  Source  locations  near  Europe  (source  locations  16-35) 
are  embedded  in  a  region  of  the  Earth  with  a  large  degree  of  structure  and  variation  in  phase 
velocity.  As  a  consequence,  the  Ms  corrections  tend  to  be  larger  than  most  other  source  locations 
on  the  map.  On  the  other  hand,  source  locations  in  central  and  northern  China  tend  to  exhibit 
smaller  Ms  corrections  because  changes  in  phase  velocity  occur  over  longer  distances,  leading  to 
smaller  gradients. 


0  20  40  60  80  100  120  140  160 


Longitude  (deg) 

Figure  52.  Map  of  the  gradient  of  the  20-second  phase  velocity  spanning  Eurasia  and  most  of  Africa 
and  Australia.  The  overlaid  scatter  plot  contains  the  simulated  explosion  source  locations,  with  size  of 
the  scatter  points  proportional  to  the  angle-averaged  rms  MS  corrections  corresponding  to  the  source 
location. 

In  addition  to  studying  the  amplitude  of  Ms  corrections,  it  is  clearly  important  to  more  fully 
understand  the  effects  of  phase  velocity  uncertainty  on  our  results.  For  example,  when  the 
uncertainty  in  the  Ms  correction  (i.e.,  standard  deviation  over  the  ensemble)  becomes  as  large  as 
the  magnitude  of  the  correction,  16%  of  the  corrections  are  of  the  wrong  sign.  In  other  words, 
when  uncertainty  is  of  the  same  magnitude  as  the  correction  itself,  applying  a  Ms  correction  will 
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actually  make  the  surface  wave  magnitude  estimate  worse  than  assuming  uniform  Earth  structure 
in  16%  of  cases. 

Figure  53  shows  a  comparison  of  the  rms  amplitude  of  Ms  corrections  to  the  uncertainty  in  the 
corrections,  averaged  over  the  88  source  locations.  Because  the  curves  intersect  at  about  2000  km, 
we  conclude  that  with  uncertainty  in  the  phase  velocity  of  2-3%,  the  validity  of  the  “dandelion” 
plots  is  out  to  a  range  of  -2000  km.  Beyond  this  range,  the  magnitude  of  the  corrections  become 
comparable  to  standard  deviation  over  an  ensemble  of  possible  true  phase  velocities. 
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Figure  53.  Blue  curve:  Root-mean-squared  magnitude  of  the  Ms  corrections  (dB),  averaged  over  all 
angles  and  source  locations,  as  a  function  of  range.  Green  curve:  Average  uncertainty  in  the  Ms 
corrections,  estimated  from  the  standard  deviation  of  the  Ms  corrections  over  an  ensemble  of  25 
trials  for  each  source  location. 
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4.7.2  Comparison  of  Data  and  Model  Variability  in  Eurasia 

The  membrane  calculation  is  an  efficient  way  to  model  surface-wave  propagation  in  a  realistically 
heterogeneous  earth  model.  Based  on  the  comparisons  with  finite  difference  calculations  it  can 
reproduce  much  of  the  complex  amplitude  behavior  arising  from  diffraction  and  multi-pathing.  As 
shown  in  the  study  of  observations  in  the  western  USA,  the  method  predicts  amplitude  variation 
on  the  order  of  what  is  observed.  However,  many  factors  contribute  to  the  observed  amplitudes 
and  it  is  difficult  to  isolate  each  one  such  that  we  can  sufficiently  improve  our  earth  model  to  better 
determine  the  source  strength. 


Figure  54  (left)  shows  the  phase  velocity  model  for  20  s  Rayleigh  waves  around  a  source  located 
in  Kyrgyzstan.  Figure  54  (right)  shows  the  log  of  ratio  of  the  maximum  amplitude  of  membrane 
waves  propagating  in  that  heterogeneous  model  to  the  maximum  amplitude  of  membrane  waves 
propagating  in  a  uniform  (homogeneous  average)  model.  Both  membrane  calculations  are  for  an 
explosion  source  to  eliminate  any  effects  from  radiation  pattern.  The  membrane  calculation 
predicts  narrow  regions  of  enhanced  and  reduced  amplitudes  with  peak-to-peak  variations  of 
nearly  a  magnitude  unit.  These  variations  occur  over  small  changes  in  azimuth. 
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Figure  54.  Left:  Map  of  20-s  Rayleigh  wave  phase  velocity  around  a  source  located  in  Kyrgyzstan. 
The  model  is  the  Leidos  model,  Stevens  et  al.  (2008).  Right:  log  of  the  amplitude  ratio  of  the 
maximum  amplitude  of  the  membrane  waves  computed  for  an  explosion  source  and  the  model 
shown  in  (a)  to  the  maximum  amplitude  of  the  membrane  waves  computed  for  an  averaged, 
uniform  model.  Note  the  predicted  fine  structure  of  strong  variation  in  amplitudes  predicted  for  the 
Leidos  model. 


One  contributor  to  the  variation  in  the  calculated  amplitudes  is  the  roughness  of  the  phase-velocity 
model,  which  in  some  cases  is  caused  by  unrealistic  boundaries  between  grid  cells.  We  learned  in 
our  earlier  work  in  the  western  USA  that  smoother  models  reduce  the  small-scale  variation  in 
amplitudes,  although  the  smaller  source-receiver  distances  in  that  study  further  reduced  the  overall 
amplitude  variations.  To  assess  the  effects  of  roughness,  we  smooth  the  model  and  recalculate 
surface  waves.  In  Figure  55,  we  derive  a  smoothed  phase-velocity  model  from  the  model  in  Figure 
54  and  repeat  the  calculation.  As  we  expected,  the  smoothed  model  does  reduce  the  small-scale 
variation.  However  the  relatively  narrow  regions  of  strong  amplification  persist  and  the  amount  of 
variation  is  still  considerable. 
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Figure  55.  Top  Left:  smoothed  Leidos  model 
centered  on  a  specific  event.  Top  Right:  Surface 
wave  amplitude  variations  in  the  smoothed 
model  (compare  with  Figure  54).  Right: 
Difference  in  amplification  between  rough  and 
smooth  models.  Red  dots  are  station  locations . 
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Comparing  the  previous  figures,  it  appears  that  the  dominant  regions  of  amplification/reduction 
persist.  However,  the  details  vary  considerably.  Also  shown  in  Figure  55  is  a  comparison  of  the 
two  models,  taking  the  difference  between  the  initial  rougher  model  and  the  smoothed  model.  The 
rougher  model  clearly  shows  much  larger  variability  as  expected,  but  it  is  also  clear  that  the  regions 
of  extreme  amplitude  variation  have  also  slightly  relocated.  The  predicted  amplitude  features  from 
similar  models  predict  an  amplification  pattern  that  is  grossly  similar,  but  that  is  very  different  in 
the  details,  details  directly  related  to  our  observational  capability. 

The  kernels  used  in  the  Born  approximation  illustrate,  at  least  qualitatively  the  issues  involved. 
The  kernels  for  phase  delay  -  i.e.,  travel  time  -  show  that  a  predicted  phase  advance  at  a  station  is 
directly  correlated  with  an  increase  in  phase  velocity  along  the  path  from  the  source,  as  one  would 
expect  from  simple  geometry.  There  are  some  regions  of  negative  correlation  (i.e.  the  negative 
side  lobes  modeled  in  an  earlier  section)  in  which  the  phase  is  advanced  by  a  decrease  in  phase- 
velocity  but  these  regions  are  relatively  inconsequential.  On  the  other  hand  the  Born  amplitude 
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kernels  are  narrower  with  deep  symmetric  side  lobes  and  show  clear  sensitivity  to  the  variation,  or 
roughness,  of  the  earth  model,  perpendicular  to  the  wave  propagation  path.  The  tomographic 
procedure  is  insensitive  to  this  structure  on  the  necessary  scale,  so  it  is  not  surprising  that  earth 
models  derived  from  phase  and  group  delays,  and  even  those  incorporating  attenuation  may  not 
predict  the  observed  amplitude  variations. 


Figure  56  shows  how  magnitude  differences  (more  precisely,  the  log  spectral  amplitude  ratios 
shown  above  in  Figure  54),  vary  with  distance  for  the  initial  and  smoothed  phase-velocity  models 
for  three  different  events  in  Eurasia.  For  each  model,  the  mean  and  standard  deviation  of  the 
amplitude  ratios  are  shown  for  annuli  1000  km  wide  between  1000  and  5000  km  from  the  source. 
At  20  s,  there  is  about  a  0.1  reduction  in  standard  deviation  between  the  initial  Leidos  model  (e.g.. 
Figure  54)  and  the  smoothed  model  (e.g.  Figure  55).  There  is  some  variation  with  source  region 
and  it  is  more  pronounced  with  the  initial  Leidos  model. 


Figure  56.  Log  amplitude  ratio  with  numerator  the  amplitude  of  an  explosion  source  propagated 
using  the  membrane  calculation  with  velocity  taken  from  the  Leidos  Eurasian  earth  model  and  the 
denominator  the  amplitude  of  an  explosion  source  propagated  on  a  uniform  membrane.  The  black 
symbols  are  for  the  smoothed  Leidos  model;  the  red  symbols  are  for  the  unaltered  Leidos  model.  The 
points  represent  averages  in  1000  km  annuli  from  1000  to  5000  km  and  surrounding  each  of  three 
source  locations.  Left:  20  s  period.  Right:  10  s  period. 

As  expected,  the  smoothed  model  reduces  the  variance  in  the  magnitude  measures  compared  to 
those  from  the  initial  model.  It  is  then  of  interest  to  see  how  well  the  variance  compares  with  that 
from  actual  magnitude  measurements. 

As  described  earlier  in  the  report  we  have  measured  amplitudes  for  a  large  number  of  Eurasian 
events  using  data  carefully  selected  for  quality.  Of  the  initial  set  of  over  1600  events,  500  were 
selected  to  have  enough  observational  data  to  make  a  reasonable  estimate  of  network  magnitude. 
We  consider  three  cases:  (1)  path  equalized  amplitudes  using  path-dependent  attenuation  resulting 
from  the  inversion  described  above;  (2)  Ms(VMAX)  magnitudes  using  the  amplitudes  from  (1)  and 
corrected  according  to  the  formula  from  Russell  (2006),  and  (3)  Ms(VMAX)  magnitudes  measured 
from  path-equalized  synthetics,  which  include  the  effects  of  the  radiation  pattern  of  the  event.  In 
each  case  we  measure  amplitudes  in  the  time-domain  using  narrow-band  Butterworth  filters 
according  to  the  Ms(VMAX)  procedure  from  Russell.  Despite  the  quality-control  screening  we 
applied,  there  remain  observations  that  appear  to  be  large  outliers,  much  larger  than  anything 
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predicted  above.  To  reduce  the  influence  of  these  measurements  we  calculate  the  network 
magnitude  for  each  event  using  the  median.  The  standard  deviation  is  then  approximated  by  the 
median  absolute  deviation,  or  MAD,  between  each  observation  and  its  associated  network 
magnitude.  Figure  57  shows  how  the  magnitude  measures  vary  with  distance  with  the 
measurements  partitioned  as  described  in  Figure  56.  The  MAD  is  corrected  to  be  a  proxy  for  the 
standard  deviation,  so  the  vertical  axes  are  directly  comparable.  The  synthetic  case  shows  how  the 
source  radiation  pattern  affects  the  variation  of  the  magnitude  measurements.  The  radiation  pattern 
is  significant  but  is  not  obviously  the  dominant  contribution.  Note  also  for  this  data  set,  the 
corrections  used  to  compute  Ms(VMAX)  appear  to  under  correct  the  magnitude  at  distance  ranges 
less  than  2000km.  The  higher  attenuation  derived  from  the  Q  inversion  discussed  earlier  provides 
a  better  distance  correction. 
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Figure  57.  Variation  of  magnitude  measurements  relative  to  their  associated 
network  averaged  magnitude.  Distance  bins  are  1000  km  wide  from  1000  to  5000  km. 
Black  symbols  are  path  equalized ,  attenuation  corrected  amplitudes.  Red  are 
magnitudes  corrected  using  the  correction  formula  in  Russell ,  2006.  Blue  symbols  are 
for  path  equalized  synthetic  data ,  which  removes  everything  but  the  radiation  pattern. 
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5.  Conclusions 


The  main  tasks  of  this  project  are: 

1.  Development  and  testing  of  a  membrane  surface  wave  code  (SurfMembrane) 

2.  Collection  of  US  Array  data  and  earth  models  in  the  area  covered  by  the  US  Array  to  test 
the  membrane  surface  wave  technique  in  an  area  with  dense  coverage  and  well-developed 
earth  models. 

3.  Collection  of  Eurasian  and  North  African  surface  wave  data  and  earth  models. 

4.  Inversion  of  Eurasian  and  North  African  surface  wave  data  for  Q  and  development  of  path 
dependent  attenuation  maps. 

All  of  these  tasks  have  been  completed.  SurfMembrane  was  developed  during  year  1,  and  during 
year  2  it  was  improved  and  the  capability  to  calculate  finite  frequency  kernels  was  added. 
SurfMembrane  3.0  is  being  delivered  together  with  this  report.  The  data  sets  of  both  US  Array  data 
and  Eurasian/ African  data  were  retrieved  from  IRIS.  Centroid  Moment  Tensor  (CMT)  solutions 
were  found  for  as  many  events  as  possible  during  the  time  period  of  interest,  and  locations  and 
depths  were  updated  with  ISC  solutions  for  the  same  events. 

We  inverted  a  large  Eurasian  and  northern  African  data  set  for  Q  structure.  The  data  set  retrieved 
from  IRIS  consisted  of  59,000  waveforms  from  1850  Eurasian  and  African  earthquakes,  all  with 
CMT  solutions.  These  were  added  to  the  smaller  data  set  used  in  the  Stevens  et  al  (2008)  study. 
Both  the  new  and  the  old  data  sets  were  subjected  to  an  intensive  quality  control  procedure  to 
remove  bad  or  questionable  data,  deep  events  and  events  with  inaccurate  CMT  solutions.  The  final 
data  set  contained  23,148  waveforms  from  998  events.  Inversion  results  are  similar  to  the  2008 
study  and  cover  a  larger  region  of  Eurasia.  There  is  a  band  of  high  attenuation  stretching  across 
the  Middle  East  from  the  Mediterranean  Sea  to  India. 

To  validate  SurfMembrane,  calculations  were  compared  with  the  results  of  a  3D  finite  difference 
calculation  corresponding  to  the  complex  region  between  Lop  Nor  and  KNET  (Stevens  et  al, 
2008).  SurfMembrane  was  run  for  each  frequency  and  compared  with  results  from  the  3D 
calculation  filtered  in  approximately  the  same  frequency  band.  The  SurfMembrane  corrections 
significantly  improved  the  consistency  of  the  results,  which  shows  that  if  the  source  and  structure 
are  well-defined,  SurfMembrane  can  be  used  to  generate  structural  amplitude  corrections. 

We  ran  SurfMembrane  for  73  events  with  CMT  solutions  within  the  US  Array  and  compared 
calculated  and  observed  amplitude  variations.  This  was  less  successful  for  two  reasons:  1)  the 
structure  is  not  known  exactly;  2)  radiation  pattern  effects  were  strong  enough  to  dominate  over 
the  structural  variations.  Because  radiation  pattern  effects  are  so  strong,  small  errors  in  the  CMT 
solution  swamp  the  variations  due  to  the  structural  effects  that  we  are  trying  to  model.  Although 
the  SurfMembrane  corrections  derived  with  an  explosion  source  had  a  small  positive  effect,  the 
improvement  is  not  very  significant. 

We  applied  SurfMembrane  to  events  in  the  Eurasian/North  African  region.  Although  the  data 
density  is  not  as  high  as  for  the  US  Array,  the  paths  are  much  longer  and  so  both  structural  effects 
and  anelastic  attenuation  are  considerably  larger.  We  ran  SurfMembrane  calculations  for  88  well- 
distributed  earthquakes  from  the  Eurasian  data  set  and  compared  the  results  with  data.  We  find 
that  the  SurfMembrane  calculations  can  be  used  to  correct  Ms  for  diffraction  out  to  distances  of 
about  2000  km,  however  at  greater  distances  diffraction  is  too  sensitive  to  small  variations  in  earth 
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structure  to  allow  reliable  predictions  for  typical  variations  and  uncertainties  in  phase  velocity 
maps. 

A  large  contributor  to  the  observed  amplitude  variation  in  the  surface  waves  is  the  ‘roughness’  in 
the  spatial  variation  of  the  earth  model.  The  scatter  in  the  observed  amplitudes,  or  equivalently, 
magnitude  corrections,  can,  therefore,  help  us  specify  smoothness  constraints  on  our  earth  models, 
which  is  an  input  to  the  tomographic  inversion.  The  overall  scatter  we  observe  in  the  magnitudes 
is  consistent  with  predictions  based  on  the  spatial  variation  seen  in  the  Leidos  1 -degree  earth  model 
and  smoother  models  predict  much  smaller  amplitude  scatter,  especially  at  the  nearer  distances. 
But  there  remains  substantial  uncertainty  in  how  radiation  pattern  is  affecting  the  scatter,  even  with 
CMT  solutions.  Careful  study  of  multiple  nearby  events  using  joint  inversion  methods  would  help 
further  isolate  the  effect  of  propagation  and  provide  a  more  realistic  estimate  of  the  scatter  in 
amplitudes  due  to  propagation  effects.  For  example,  the  relative  homogeneity  of  the  Leidos  model 
for  large  areas  of  Asia  as  manifest  in  Figure  52  predicts  that  scatter  in  amplitudes  due  to 
propagation  should  be  relatively  small  in  this  region.  Further  study  including  the  effects  of  the 
source  could  verify  this  prediction.  The  SurfMembrane  approach  provides  an  efficient  way  to 
incorporate  amplitude  constraints  into  inversion  studies  and  make  realistic  predictions  from  the 
results. 
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Appendix  A.  SurfMembrane  Description  and  User  Manual 
A.l  Input  to  SurfMembrane 

SurfMembrane  is  written  in  Fortran  95  with  some  Fortran  2003  extensions.  Input  to  the 
membrane  surface  wave  code  is  in  Fortran  namelist  format  using  a  namelist  “MEM_IN”. 
The  code  is  executed  with  the  command: 

SurfMembrane  NamelistFilename 

where  NamelistFilename  is  the  name  of  a  file  containing  the  input  namelist. 

An  example  input  file  is  given  below  for  the  Lop  Nor  to  KNET  calculation  at  10  seconds 
period. 

&MEM_IN 
nx  =  678 
ny  =  858 
dx  =  2000. 
nt  =  5000 
dt  -  0.20 
printcycle  =  5 
idrive  =  302 
jdrive  =  750 
drive_type  =  2 

drivef ilename  =  ’ drive/evl . 10 . drv ’ 
cvelfilename  =  ’ cvel/KNET . 10' 
out filename  =  ’ out file /KNET . 10 . dat ’ 
stalocf ilename  =  ’ staloc/stalist ’ 
stadataf ilename  =  ’ outf ile/evl . 10 . stadata 1 
/ 

Definitions: 

nx  =  number  of  points  in  the  X  direction 
ny  =  number  of  points  in  the  Y  direction 

dx  =  spacing  between  points  (X  and  Y  spacing  is  the  same) 
nt  =  number  of  time  steps  to  run 

dt  =  time  step  (Courant  condition  applies,  must  be  <0.7  * 
minimum  travel  time  across  a  zone.  0.5  recommended) 
printcycle  =  increment  between  print  cycles.  5  will  print  cycles 
5,  10,  15  ... 

idrive  =  X  node  location  of  drive 
jdrive  =  Y  node  location  of  drive 
drive_type  =  0  No  drive 

1  Not  currently  implemented 

2  Specified  fixed  velocity 

3  Specified  velocity  added  to  calculated  velocity 
drivef ilename  =  drive  file  name.  See  Appendix  for  format 
cvelfilename  =  phase  velocity  map  file  name.  See  Appendix  for 

format 

outf ilename  =  output  filename.  See  Appendix  for  format 
stalocf ilename  =  station  location  filename  (saved  every  cycle) . 

See  Appendix  for  format 
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stadataf ilename  =  output  station  time  series.  See  Appendix  for 
format . 


Other  namelist  parameters  not  used  in  the  example  above,  with  default  values: 


cvel_in  =  1.0.  Constant  phase  velocity  value  used  if  cvelfilename 
is  not  input. 

initstatef ilename  =  none.  Specify  initial  wavefield  in  grid. 
force_type  =  0 .  No  forcing  function 

=  1.  Calculated  Gaussian  forcing  function 
=  2 .  Specified  forcing  function 

forcef ilename  =  none.  Forcing  function  filename.  See  Appendix  for 
format . 

printspacing  =  1.  Spacing  between  points  on  output  file  (3  means 
every  third  point  is  printed) . 

reverse_f orce  =0.  If  positive,  time  reverse  the  input  force 
(used  for  adjoint  calculations) . 
isrc  =  0.  X  location  of  forcing  function 
jsrc  =0.  Y  location  of  forcing  function 
Forcing  function  variables  (see  below) 
tau  =  20.0 

tauO  =  2.628 

tstart  =  48 . 

input_only  =0.  If  positive,  print  out  input,  and  then  stop 
without  running  the  calculation. 


The  calculated  forcing  function  is  illustrated  in  Figure  A-l 
anddefined  by: 


a  = 


A  — 


A 

T 

-2a3 


t  =  time  -  tstart 

F(t)  =  At2  expf-(atf)2  j 


Figure  A-l.  Default  forcing  function. 
Input  variables  for  calculation  of  sensitivity  kernels: 


kernel_f rom_stadata  =0.  If  positive,  calculate  kernel  from 
stadata  file. 

kernelf ilename  =  'none' .  Output  kernel  file. 

sourcememf ilename  =  'none' .  Input  file,  results  of  previous 

membrane  calculation.  Not  needed  if  kernel_f rom_stadata>0 . 
In  that  case  outf ilename  is  used  as  input  for  kernel  calc. 
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A.2  SurfMembrane  File  Formats 


In  this  appendix  we  describe  the  formats  needed  for  reading  and  writing  the  files  used  by 
SurfMembrane.  Error  checking  in  the  open  and  read/write  statements  has  been  removed 


for  clarity.  In  the  examples  in  this  appendix,  the  X  axis  has  been  chosen  to  be  north  and  the 


Y  axis  east,  because  these  were  the  coordinates  used  in  the  earlier  3D  KNET  calculation. 
This  is  not  a  requirement  and  X  and  Y  can  be  chosen  to  be  east  and  north  as  long  as 
coordinates  are  consistent  in  all  input  files.  Note  that  data  on  binary  files  are  stored  in 
column  order,  so  all  rows  are  read  sequentially  by  column  for  each  matrix. 

A.2.1  Phase  Velocity  File  (input)  cvelfilename.  Fortran  direct  access 
binary  file. 

Opened  with: 

open(VELFILE,file=cvelfilename,form=,unformatted’,access=,direct’,status=,old,,recl=nx 

*ny*4) 

Read  with: 

read(VELFILE,rec=l)  cvel(l  :nx,l  :ny) 

The  image  in  Figure  A-2  was  made  with  the  following  Matlab  commands: 

FILE=f open ( cvelfilename , ’ r ’ ) ; 
cvel=f read (FILE, [nx  ny] , ’ f loat32 ’ ) ; 
f close  (FILE) ; 
image sc (y , x, cvel ) ; 

axis  equal;  set (gca, ’ Ydir Normal ’) ; 
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Figure  A-2. 10  second  phase  velocity  map  used  in  KNET  calculation. 
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A.2.2  Output  File  outfilename  (output).  Fortran  direct  access  binary 
file. 

Opened  with: 

open(OUTFILE,file=outfilename,form='unformatted',access='direct',status='replace', 
recl=(  (nx+2)*(ny+2)*4) ) 

Note  that  one  additional  row  and  column  are  added  to  each  of  the  sides,  top  and  bottom, 
so  the  dimensions  of  the  output  file  are  nx+2  and  ny+2  instead  of  nx,ny 

Written  with  (in  time  loop): 

if(  mod(it,printcycle)  ==  0)  then 

write(OUTFILE,rec=it/printcycle)  ut(0:nx+l,0:ny+l) 
endif 

The  image  in  Figure  A-3  was  made  with  the  following  Matlab  commands,  where  nt  was 
chosen  to  be  the  number  of  saved  states  corresponding  to  240  seconds. 

xx=0 : 2 : nx*2 ; 
yy=0 : 2 : ny*2 ; 
for  i=l : nt 

Z  =  f read (FILE ,  [nx+2  ny+2 ] ,  ’ float 32 ’)  ; 

end 

xs=302 ; ys=750 ;  %Source  location 

Z (xs-4 : xs+6 , ys-4 : ys+6 ) =0 ;  %Blank  out  source  region 

image  s  c ( y y , xx , Z ) ; 

axis  equal; colorbar; 

set (gca,  ’ Ydir ’ ,  ' Normal ’ ) ; 


Figure  A-3.  Snapshot  of  velocity  from  KNET  calculation  at  240  seconds. 
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A.2.3  Station  Location  File  (input)  stalocfilename.  Ascii  file. 

Opened  with: 

open(STALOCFILE,file=stalocfilename,form='formatted',status='old') 
Read  with: 

read(STALOCFILE,*,iostat=ios,iomsg=message)  nstation 
do  ista=l, nstation 

read(STALOCFILE,*)  istation(ista),jstation(ista) 
end 


Here  is  an  example  of  a  station  location  file  from  the  KNET  run: 
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A.2.4  Station  Data  File  (output)  stadatafilename.  Ascii  file. 

Opened  with: 

open(STADATAFILE,file=stadatafilename,form='formatted',status='replace') 
Written  with  (in  time  loop,  nstation  columns  in  esl5.5  format): 
write(formatstring,'(al,i2,a)'),'(',nstation,'esl3.5)' 

write(STAD AT AFILE, formatstring)  (ut(istation(ista)j  station(ista)),ista=  1 , nstation) 
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A.2.5  Initial  State  File  (input)  initstatefilename.  Fortran  direct  access 
binary  file. 

Opened  with: 

open(INITSTATEFILE,file=initstatefilename,form='unformatted',access='direct',status 

old',recl=nx*ny*4) 

Read  with: 

read(INFILE,rec=l)  ut(l:nx,l:ny) 
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A.2.6  Drive  File  (input)  drivefilename.  Ascii  file. 

Opened  with: 

open(DFILE,file=filename,form='formatted',status='old') 
Read  with: 

read(DFILE,*)  nxydrive,ntdrive 
read(DFILE,  *)  (tdrive(itdrive),itdrive=  1  ,ntdrive) 
do  ixydrive=l,nxydrive 

read(DFILE,*)  ixdrive(ixydrive),iydrive(ixydrive) 
read(DFILE,*)  (drive(itdrive,ixydrive),itdrive=l  ,ntdrive) 
end  do 
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A.2.7  Forcing  File  (input)  forcefilename.  Ascii  file 

Opened  with: 

open(FFILE,file=filename,form='formatted',status='old') 

Read  with: 

read(FFILE,  *  ,iostat=ios,iomsg=message)  nforce 
do  if orce=l, nforce 

read(FFILE,*,iostat=ios,iomsg=message)  tforce(iforce),fforce(iforce) 
end  do 


Approved  for  public  release;  distribution  is  unlimited. 
66 


A.2.8  Sensitivity  Kernel  (output)  File  Format 

Opened  with: 

open(KRNFILE,file=kernelfilename,form= 'unformatted', access='direct',status='replace', 
recl=(nx+2)*(ny+2)*4,  iostat=ios,iomsg=message) 

Written  with: 

do  ista  =  l,nstation 

write(KRNFILE,rec=ista,iostat=ios,iomsg=message)  kernelFinal(0:nx+l,0:ny+l) 

end  do 
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Appendix  B.  SurfMembrane  Code  and  Data  Deliverable 
B.l  Data  Deliverable 

The  source  code,  sample  results  and  test  data  are  all  included  in  the  file 
SurfMembrane_2017-05.tgz  which  is  distributed  with  this  report.  To  run  the  test  case,  untar 
the  file  into  the  desired  installation  directory  and  follow  the  instructions  in  the  README 
file.  The  procedure  is  summarized  below: 

1 .  Modify  the  global  Makefiles  to  have  the  correct  paths  to  compilers  installed  on  the 
system  (see  README  for  details). 

2.  Run  “make”  to  compile  SurfMembrane.  The  resulting  executable  will  be  src/ 
SurfMembrane. 

3.  Run  “make  test”  to  run  the  test  cases. 

4.  To  make  plots,  open  Matlab  from  the  installation  directory,  <distdir>,  and  enter: 
addpath  <distdir>/matlab 

addpath  <distdir>/test/matlab 
cd  <distdir>/test 
plot_examples 

5.  This  will  make  a  set  of  png  files  in  directory  test/plots  that  should  replicate  one-for- 
one  plots  in  test/refplots.  There  may  be  small  differences  due  to  differences  in  fonts 
on  different  systems. 

6.  You  can  also  run  movies  of  the  calculation  using 
“plot_example_membrane_movies”. 
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B.2  MATLAB  Scripts 

Several  MATLAB  functions  are  provided  to  simplify  data  access  and  to  make  plots  of 
calculation  results.  These  are  summarized  below.  Do  a  “help”  on  the  matlab  directory  after 
installation  for  details. 

get_max_membrane:  Retrieve  maximum  amplitudes  and  amplification  for  heterogeneous  and 
homogenous  membrane  calculations. 

[seis_max,seisuni_max,seisamp]=  get_max_membrane(memfile,memfileuni,nx,ny,nt,dt,fc,delf) 
get_membrane:  Retrieve  3D  (x,y,time)  array  of  data  from  membrane  calculation. 

memdata=get_membrane(memfile,nx,ny,nt,tlim) 

get_stadata:  Retrieve  all  calculated  station  data  from  a  membrane  calculation. 

[stadata,nsta,nt]=get_stadata(stadatafile,dt,fc,delf) 

plot_drive:  Plot  source  drive  functions. 

[h,hax]  =  plot_drive(drivefilename,plotindividual) 
plot_membrane:  Plot  images  of  spatial  data  retrieved  by  get_membrane  or  get_max_membrane. 

[h,hax]  =  plot_membrane(memdata,dxy,clim,tpause,thm,titletext) 

plot_stadata:  Plot  calculated  time  series  for  an  individual  station. 

[h,hax]  =  plot_stadata(stadata,ista,dt) 
plot_phasevel:  Plot  phase  velocity  map 
[h,hax]  =  plot_phasevel(velfile,nx,ny,dxy,xyorigin,axlim,clim) 
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B.3  Test  Case 


A  Test  Case  is  provided  that  includes  calculations  at  10  and  20  seconds  for  an  explosion  at 
0.5  km  depth  in  two  earth  structures,  corresponding  to  a  source  location  in  northeastern 
Nevada  at  latitude  41.150  and  longitude  -114.870.  This  is  the  location  of  one  of  the 
earthquakes  (event  17)  in  the  US  Array  study,  and  the  test  case  is  for  an  explosion  at  the 
same  coordinates  with  receivers  at  the  same  US  Array  locations.  Several  Matlab  scripts  are 
provided  to  make  plots  from  the  test  case.  See  section  B.l  for  instructions  on  running  the 
test  case. 


The  Matlab  script  “plot_example_phasevel”  plots  the  phase  velocity  maps  used  in  the 
calculation.  These  are  shown  in  Figure  B-l. 
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Figure  B-l.  Phase  velocity  maps  for  test  case.  Top  row  Leidos  structures ,  bottom  row 
WUS  structures .  Left  column  10  seconds,  right  column  20  seconds . 
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The  Matlab  script  “plot_example_drives”  plots  the  applied  source  functions.  These  are 
shown  in  Figure  B-2.  Each  figure  has  40  curves  and  although  they  appear  to  be  nearly 
the  same,  the  small  amount  of  phase  and  amplitude  difference  between  the  curves  is 
important.  Also,  notice  that  the  source  function  amplitudes  of  the  two  structures 
are  different,  corresponding  to  the  difference  in  surface  wave  excitation  function  for 
each  structure.  These  functions  correspond  to  the  vertical  velocity  of  the  fundamental 
mode  surface  wave  calculated  on  a  closed  surface  around  the  origin. 


Figure  B-2.  Drive  functions.  Top  row  Leidos  structures,  bottom  row  WUS  structures. 
Left  column  10  seconds,  right  column  20  seconds. 
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The  script  “plot_example_membrane”  plots  the  propagation  on  the  grid.  The  script  first 
runs  a  movie  of  each  calculation  and  then  makes  at  plot  at  150  seconds.  These  are 
shown  in  Figure  B-3. 


Figure  B-3.  Spatial  plot  of  velocity  at  150  seconds.  Top  row  Leidos  structures,  middle  row 
WUS  structures,  bottom  row  uniform  structures.  Left  column:  10  seconds,  right  column  20 
seconds. 
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The  script  plot_example_stations  plots  waveforms  at  selected  stations.  Note  that  while  the 
full  grid  is  only  saved  every  10  time  steps  (2.0  second  intervals),  the  station  data  is  saved 
for  every  time  step  (0.2  seconds),  so  the  station  data  has  higher  time  resolution.  The  plots 
made  by  this  script  are  shown  in  Figure  B-4. 


Figure  B-4.  Waveforms  at  station  #420  (location  854, 271  in  the  grid).  Top  row:  Leidos 
structures,  middle  row :  WUS  structures,  bottom  row:  uniform  structures.  Left:  10  seconds, 
right:  20  seconds. 
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The  Matlab  script  “plot_example_membrane_amp”  plots  the  spatial  distribution  of 
amplification  (or  deamplification)  caused  by  propagation  through  the  earth  structures. 
The  plots  made  by  this  script  are  shown  in  Figure  B-5. 


Figure  B-5.  Surface  wave  amplification  for  heterogeneous  relative  to  uniform  structures. 
Top  row:  Leidos  structures,  bottom  row  WUS  structures.  Left:  10  seconds,  right:  20  seconds. 
Amplification  is  in  dB. 
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Appendix  C.  Surface  Wave  Q  Inversion  Data  Deliverable 

This  final  report  is  accompanied  by  a  data  deliverable  consisting  of  the  following  files: 

LP_2017_May.tar.gz:  This  is  the  complete  set  of  global  earth,  dispersion  and  attenuation 
models.  The  contents  are  described  in  the  file  “README.models.”  It  also  contains  a 
program  compiled  under  Linux  to  calculate  the  dispersion  and  attenuation  between  any  two 
points  at  an  input  frequency. 

gamma_data_new:  These  are  all  of  the  attenuation  measurements  used  by  Leidos  during 
this  project.  The  format  of  the  file  is  given  in  the  file  “dataformat”. 

moments.txt:  This  file  gives  the  adjustments  to  CMT  moments  found  during  the  Q 
inversion  process  as  described  in  Section  4.6  of  the  report.  The  columns  correspond  to 
event  number,  natural  log  of  the  correction,  and  value  of  the  correction. 

The  data  given  in  the  datafile  are  the  original  measurements  derived  using  CMT  moments 
and  should  be  adjusted  using  the  moment  corrections  by: 

Ay  =  -AlnM0  /  r 

where  r  is  the  source  to  receiver  distance. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


AFRL 

Air  Force  Research  Laboratory 

AFSPC 

Air  Force  Space  Command 

CMT 

Centroid  Moment  Tensor 

GSN 

Global  Seismic  Network 

IRIS 

Incorporated  Research  Institutions  for  Seismology 

KNET 

Kyrgyz  Seismic  Network 

wus 

Western  United  States 
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